3

I have a problem which I try to solve with mathematica. I am having a list with x and y coordinates from a position measurement (and also with z values of the quantity which was measured at each point). So, my list starts with list={{-762.369,109.998,0.915951},{-772.412,109.993,0.923894},{-777.39, 109.998, 0.918108},...} (x,y,z). Out of some reasons, I have to fill all these x,y, and z-values into a matrix. That would be easy if I have for each y-coordinate the same amount of x-coordinates (lets say 80), then I could use Partition[list,80] which produces a matrix with 80 columns (and some rows whose number is given by the number of y-coordinates with the same value).
Unfortunately, it is not so easy, the number of x-coordinates for each y is not strictly constant, as can be seen from the attached ListPlot. xy coordinates with ListPlot Can anybody give me some suggestions, how I could fill each point of this plot / each x-y-(and z-) coordinate of my list into a matrix?

To explain better what I want to have, I indicated in the attached picture a matrix. There one can see that almost every point of my plot would fall into a cell of a matrix, only some cells would stay empty. I used in the plot the color red for the points whose x coordinates are ascending in my list and blue for the points whose x coordinate are descending in my list (the positions are measured along a meander line). Perhaps this kind of order can be useful to solve to problem... Here a link to my coordinates, perhaps this helps.

Well, I hope I explained my question well enough. I would appreciate every help much!

partial81
  • 495
  • 6
  • 14
  • It'd be nice if you could provide the dataset. If I got it, you need the matrix elements to be `z` values, and the matrix indices should approximately correspond to `x` and `y`, right? (Also, it's better not to use variable names starting with capital letter to avoid conflict with builtins/packages. `List` is a builtin!) – Szabolcs Jul 26 '11 at 15:56
  • Thanks for the quick reply. Sorry for the capital letter in list, I used it to make it better readable but did not think at the conflict problem (I edited). Well, I would love to provide the dataset here but I do not know how to upload a file... Sorry! Therefore, I hope it is ok, when I am sending a file to your email address. Yes, the matrix indices should approximately correspond to x and y. But my aim is to have (x,y,z) in each matrix element, I do not want to use interpolation of my values. – partial81 Jul 26 '11 at 16:43
  • 1
    What do you want to do with your matrix? – Dr. belisarius Jul 26 '11 at 16:47
  • @partial81 You can upload the datafile e.g. to http://ge.tt/ and include a link – Szabolcs Jul 26 '11 at 18:48
  • @Szabolcs Thanks a lot for the email reply and your answer here! Your suggestion is great! After some hours of sleep I will try to understand it better ;-) Thanks also for the hint with the upload. I did not know this page. I edited my question again and provide the link to a file with my x- and y-coordinates now. So other user can use them in the case they want to try applying yours or others suggestion. – partial81 Jul 26 '11 at 21:25
  • @belisarius First, I want to find a general way to solve my problem (just filling the coordinates into a matrix). Secondly, if I have in every cell of my matrix/tensor the x,y,z values, then I can plot my z-values with ReliefImage or MatrixPlot and not anymore with ListDensityPlot (as much as I know, with ListDensityPlot my z-values are interpolated on a grid. This is something I do not want to have all the time). Thirdly, I am sometimes interested in special regions of interest in my data/plot. Thanks to the great help of forum users, I can easily extract them if I can produce a ReliefImage. – partial81 Jul 26 '11 at 21:50
  • @partial81 If your goal was just getting a grid of `z` values, this was an awful lot of work for that... Why don't you use `Interpolation` on the `{x,y,z}` data then use the resulting interpolation function to generate a matrix? You'll be able to plot that with `ReliefImage` or whichever similar function you like. Is there any good reason why any kind of interpolation is out of question? Since the lattice is not perfectly square, my answer below also introduces distortions. – Szabolcs Jul 27 '11 at 11:37
  • @Szabolcs Ok, that is a good question. Well, first I was interested if there is a solution in general. Secondly, as I denoted in my email, I tried to find a way to extract my real values in special regions of interest out of my list with an easy method (because out of some physical reasons I sometimes may only use them). Therefore, I still think/hope that my question is useful, although it was awful lot of work for you. Perhaps others can profit by our answer one day too and appreciate your work as much as I do. – partial81 Jul 28 '11 at 08:32
  • @Szabolcs I have to admit that I tried to interpolate my slightly irregular data with mathematica, but I failed there and cannot find useful info on the web. The odd thing is, I can plot my data for example with ListPlot3D[data,InterpolationOrder->1] and get surfaces which are interpolating my data well. Now I would be interested in the interpolated values which are not exactly at the measured points (to plot e.g. cross sections). Do you know how to do this without spending much work again (or a webpage with good information)? Or do you think that this question is good enough for this forum? – partial81 Jul 28 '11 at 10:57

1 Answers1

6

The basic idea behind this solution is:

  • all points seem to lie on a lattice, but it's not precisely a square lattice (it's slanted)
  • so let's find the basis vectors of the lattice, then all (most?) points will be approximate integer linear combinations of the basis vectors
  • the integer "coordinates" of the points along the basis vectors will be the matrix indices for the OP's matrix

(The OP emailed me the datafile. It consists of {x,y} point coordinates.)

Read in the data:

data = Import["xy.txt", "Table"];

Find the nearest 4 points to each point, and notice that they lie about distance 5 away both horizontally and vertically:

nf = Nearest[data];

In:= # - data[[100]] & /@ nf[data[[100]], 5]

Out= {{0., 0.}, {-4.995, 0.}, {5.003, 0.001}, {-0.021, 5.003}, {0.204, -4.999}}

ListPlot[nf[data[[100]], 5], PlotStyle -> Red, 
  PlotMarkers -> Automatic, AspectRatio -> Automatic]

enter image description here

Generate the difference vectors between close points and keep only those that are about length 5:

vv = Select[
      Join @@ Table[(# - data[[k]] & /@ nf[data[[k]], 5]), {k, 1, Length[data]}], 
      4.9 < Norm[#] < 5.1 &
     ];

Average the vectors out by directions they can point to, and keep two "good" ones (pointing "up" or to the "right").

In:= Mean /@ GatherBy[vv, Round[ArcTan @@ #, 0.25] &]

Out= {{0.0701994, -4.99814}, {-5.00094, 0.000923234}, {5.00061, -4.51807*10^-6},  
      {-4.99907, -0.004153}, {-0.0667469, 4.9983}, {-0.29147, 4.98216}}

In:= {u1, u2} = %[[{3, 5}]]

Out= {{5.00061, -4.51807*10^-6}, {-0.0667469, 4.9983}}

Use one random point as the point of origin, so the coordinates along the basis vectors u1 and u2 will be integers:

translatedData = data[[100]] - # & /@ data;

Let's find the integer coordinates and see how good they are (how far they are from actual integers):

In:= integerIndices = LinearSolve[Transpose[{u1, u2}], #] & /@ translatedData ;

In:= Max[Abs[integerIndices - Round[integerIndices]]]

Out= 0.104237

In:= ListPlot[{integerIndices, Round[integerIndices]}, PlotStyle -> {Black, Red}]

enter image description here

All points lie close to the integer approximations.

Offset the integer coordinates so they're all positive and can be used as matrix indices, then gather the elements into a matrix. I put the coordinates in a point object in order not to confuse SparseArray:

offset = Min /@ Transpose[Round[integerIndices]]
offset = {1, 1} - offset

result = 
 SparseArray[
  Thread[(# + offset & /@ Round[integerIndices]) -> point @@@ data]]

result = Normal[result] /. {point -> List, 0 -> Null}

And we finally have a matrix result where each element is a coordinate-pair! (I was sloppy doing 0 -> Null here to mark missing elements: it's important that data contained no exact 0s.)

MatrixForm[result[[1 ;; 10, 1 ;; 5]]]

enter image description here

EDIT

Just for fun, let's look at the deviations of points from the precise integer lattice sites:

lattice = #1 u1 + #2 u2 & @@@ Round[integerIndices];

delta = translatedData - lattice;
delta = # - Mean[delta] & /@ delta;

ListVectorPlot[Transpose[{lattice, delta}, {2, 1, 3}], VectorPoints -> 30]

enter image description here

Szabolcs
  • 24,728
  • 9
  • 85
  • 174
  • Nice idea with finding the basis vectors. What breaks if you do not use point with data in the sparse array? (writing from iPad, so no mathematica to try) – acl Jul 26 '11 at 19:43
  • @acl, I didn't spend time to figure out how it works exactly ... what I knew before is that one can create a (sparse) array with `SparseArray[{ {row, column} -> element }]`. If `element` is a `List`, `SparseArray` gives me a 1D array. I guess it's because `SparseArray[{index, index, ...} -> {element, element, ...}]` is also a valid syntax. – Szabolcs Jul 26 '11 at 19:46
  • @Sza Curious! I did almost the same thing here http://stackoverflow.com/questions/4917896/recognizing-distortions-in-a-regular-grid/4917954#4917954 – Dr. belisarius Jul 26 '11 at 22:21
  • @Szabolcs Thank you for your solution and the explanations. They are really fine! It is just a bit difficult for me to understand all details because your code is far above my programming skills. But anyway, I got the idea and understand much of the solution! And it seems that I can adapt your solution well to my real data. Thank you again! – partial81 Jul 27 '11 at 12:39