coordinatesqgisshapefile
I have a project in QGIS with two layers, one with several geographical zones (imported from a Shapefile file), and one with several points (imported from a CSV with longitude/latitude coordinates).
I would like to map each point to the ID of the zone that contains that point, i.e., create a new CSV in the form "longitude,latitude,ID". How can I do it in QGIS (or, possibly, any other "scripted" solution)?
You have two tasks: (a) find the grid data for the vertices of a grid cell in which a GPS location (x,y) lies and (b) interpolate those data.
An efficient way to accomplish (a) is to store the grid as an array with a single index. Often grids are stored by rows with data progressing left to right within each row. A grid is determined by the coordinates of its origin, (x0,y0) (which will be (-180, -90) for a worldwide grid in decimal degrees), its cellsize c (which is 1/24 in your case), and the number of points per row (n). The values would be recorded in the sequence
(x0,y0), (x0+c,y0), (x0+2c,y0), ..., (x0+(n-1)*c,y0)
(x0,y0+c), (x0+c,y0+c), (x0+2c,y0+c), ..., (x0+(n-1)*c,y0+c)
...
(x0,y0+(m-1)*c), (x0+c, y0+(m-1)*c), ..., (x0+(n-1)*c,y0+(m-1)*c)
This sequence can be addressed with a single index k = 0, 1, ..., n*m-1.
Given an arbitrary point (x,y), you use simple calculations to find the indexes where its four nearest neighbors are stored: they lie between rows i and i+1 and columns j and j+1 where
i = Floor((y-y0)/c), j = Floor((x-x0)/c).
The corresponding offsets are ni + j, n(i+1) + j, ni + (j+1), and n(i+1) + (j+1).
This assumes x0 <= x < x+nc and y0 <= y < y+mc, which is worth checking in a general-purpose program. The coordinates of these four points are readily seen to be equal to
(x0 + j*c, y0 + i*c), (x0 + j*c, y0 + (i+1*c), etc.
Now that you can solve (a) efficiently, it remains to interpolate the values found, which I will denote z00, z01, z10, and z11, situated as suggested in this crude picture (and shown more precisely in the figure at the end):
z01 z11
z
z00 z10
You wish to estimate the value of z. It is convenient to use bilinear interpolation. The formulas amount to the following: first interpolate linearly along the top and bottom edges of the square. Here, only the x coordinate is changing (from x0+j*c to x0+(j+1)*c). Rescale this to t varying from 0 to 1, where
t = (x - (x0+j*c))/c = (x-x0)/c - j.
The interpolated values are
z0 = (1-t)*z00 + t*z10, z1 = (1-t)*z01 + t*z11
Now interpolate vertically between z0 and z1. First rescale y to s varying from 0 to 1, where
s = (y-y0)/c - i
and find the interpolated value
z = (1-s)*z0 + s*z1.
Let (x0,y0) = (-180,-90), c = 1/24 = 0.041667, n = 24*360 = 8640, and (x,y) = (54.4786674627, 17.0470721369) (that is, near latitude 17 degrees north and longitude 54.5 degrees east). Then
i = Floor((y-y0)/c) = Floor(24*(90 + 17.047...)) = 2569
j = Floor((x-x0)/c) = Floor(24*(180 + 54.478....)) = 5627.
The undulation for (i,j), namely z00, is found at offset 8640*i+j = 48,619,849. The undulation for (i+1,j), namely z01, is found at offset 8640*(i+1)+j = 48,628,489. The undulation for (i,j+1), namely z10, is found at offset 48,619,850, and the undulation for (i+1,j+1), namely z11, is found at offset 48,628,490. Suppose these have the values given; that is,
z00 = 31.945, z01 = 31.993, z10 = 31.866, and z11 = 31.911
To interpolate, compute
t = (x-x0)/c - j = 0.488019
s = (y-y0)/c - i = 0.129731
Therefore
z0 = 31.96842 and z1 = 31.88796
Finally
z = 31.958.
From the little information you give, it is hard to tell what goes wrong.
Most common errors are:
To test what is right and what is wrong, use a background layer of the Openlayers plugin, or Natural Earth shapefiles.
Best Answer
You need to do the following:
right-click the file in the table of content -> Save as -> Choose shapefile formatand make sure you select the CRS that matches the polygon shapefile CSRVector -> Data Management Tools -> Join attribute by locationright-click the file shapefile in the table of content -> Save as -> Choose CSV file. In thelayer Option -> Geometry -> choose 'Default'if you want to keep the original coordinate without appending another XY field to CSV output file, as you can see below: