3

Initial problem: I am currently trying to create a spatial grid in which each cell is centered on predefined spatial points. These spatial points are projected in the WGS84 international coordinate system.

I am starting with a long/lat dataframe containing around 8000 point-location over France (data available here).

df=read.table("~/centroids.csv",sep=",",header=T)
plot(df$longitude,df$latitude,pch=".")

centroidsPlotFrance

I want to compute "tile" polygons that will fill the entire space around these points. Something like in this post, with each point of df at the center of tiles, such as :

exampleTileCenteredOnPoints

I tried several methods from the sf and raster package such as : st_make_grid(), rasterFromXYZ() combined with rasterToPolygons() that I saw in this post, or even the basic raster() using gridded() that I saw on this post. I also tried to use gBuffer() such as in this post. But without any sucess... I am always fighting against projection or non-regular/planar coordinates errors !

Once I will be able compute those tiles (or polygons), I want to store them in a simple feature object (from the sf package). Within this sf object: each polygon will correspond to a cell of the final spatial grid; each polygon will be centered on a spatial point from df; and each polygon will inherit the ID of the point (from df) on which it is centered.

Is that possible ???

3
  • 1
    Those points don't look to be on a regular grid - look at just Corsica for example - are you expecting a regular polygon grid as output? Those points might be on a regular grid in another coordinate system but it seems to break down in places.
    – Spacedman
    Jul 8 at 19:26
  • I agree, points are not on regular grid. But i wanted to find a way to shape polygon that could fill the space around each point (the "grid" will not be regular for sure). Do you think this is impossible ? Maybe I have to find the one projection in which those points are in a regular grid ? Maybe this projection does not even exists... Jul 8 at 20:06
  • 1
    I tried with different coordinate systems, and I found out that your points, when displayed in epsg:27572 (NTF Lambert Zone II), form a regular grid of 8km x 8km, with at most around 2-3 meters of error across the whole country. So I guess generating your grid in epsg:27572 would overlay quite well.
    – FSimardGIS
    Jul 8 at 22:33
0

Final solution:

library("sp")
library("sf")
library("raster")
library("parallel")
library("ggplot2")

df=read.table("~/centroids.csv",sep=",",header=TRUE)

df2=df
coordinates(df2)=~longitude+latitude
proj4string(df2)=CRS("+init=epsg:4326") # points are initially projected in WGS84
df2=spTransform(df2,CRS("+init=epsg:27572")) # points are distributed on an approximatly regular grid in NTF Lambert Zone II (thank you @FSimardGIS !!!)

grid=points2grid(df2,tolerance=0.00587692) # 0.00587692 is the minimum value I found by hand, by trial and error
grid=as.data.frame(SpatialGrid(grid)) 
coordinates(grid)=~longitude+latitude

grid=rasterFromXYZ(grid)
grid=rasterToPolygons(grid)
grid=st_as_sf(grid)
st_crs(grid)=27572
grid=st_transform(grid,crs="+init=epsg:4326") # re-projection of the grid in the initial projection of df values which are in WGS84

cl=makeCluster(spec=(detectCores()-1)) # parallel computation in order to know which cells from grid actually contain points from df
clusterEvalQ(cl,library("sf"))
clusterExport(cl,c("grid","df"))
df$gridLoc=parApply(cl,df,1,function(row){
  pnt=st_sfc(st_point(c(row["longitude"],row["latitude"])),crs=4326)
  id=row["ID"]
  return(as.numeric(st_intersects(pnt,grid)))
})
stopCluster(cl)

grid=grid[df$gridLoc,] # selection of cells from grid that actually contain points from df
grid$layer=df$ID # transfert the ID of centroids to the coresponding cell of the grid

x11()
ggplot()+
  geom_sf(aes(fill=layer),col="black",data=grid)+
  geom_point(aes(x=longitude,y=latitude),pch=21,fill="white",col="black",size=.9,data=df)+
  scale_fill_gradient2(low="#47fe44",mid="#f8fe44",high="#fe4444",midpoint=5000)

Which give this nice plot (color layer corresponds to the initial ID of spatial points in df, white dots correspond to the initial centroids from df): enter image description here

Even if some points seem to be a bit off-center from the cell, its only a plotting artifact due to the high resolution, as we can see here with increased zoom on Corsica: enter image description here

The answer to this post is therefore : (1) find the right projection in which points form a regular grid; (2) use gridded() or points2grid() to compute the regular grid; (3) transform the grid points into a raster using rasterFromXYZ(); (4) transform the raster into polygons using rasterToPolygons(); (5) re-project the polygons into the initial projection using st_transform(); (6 optional) if you do not want to keep all created cells within the grid, use st_intersects() to find which cells actually contains initial points.

Thanks a lot for your advices !

Testing the grid centroid offsets: I can compute the centroids of the grid and compare their location to the initial centroids:

estimatedCentroids=st_centroid(grid) # new centroids from grid

coordinates(df)=~longitude+latitude
proj4string(df)=CRS("+init=epsg:4326")
trueCentroids=st_centroid(st_as_sf(df)) # initial centroids from df

centroidsOffset=st_distance(estimatedCentroids,trueCentroids,by_element=TRUE) # compute distance of each pair of centroids between grid and df

hist(centroidsOffset,nclass=100)

With the final solution that I present above, I therefore make an average error of 550m (range from 170m to 900m): enter image description here

If I well understand the methodology, this is due to the use of the points2grid() function with a tolerance parameter of 0.00587692. Which means that points are not distributed on a true regular grid, even in the NTF Lambert Zone II projection (which still the better I found).

4
  • 1
    Interesting. When I tested it on my side, the errors were more in the range of 2-3 meters. The points really closely matched a regular grid. But I generated my grid from scratch in epsg:27572, not with the points to grid function. The grid was 8000 x 8000, with a range of (0, 1240000, 1605000, 2685000)
    – FSimardGIS
    Jul 9 at 17:46
  • Could this difference be due to the fact that I computed distances in the initial projection (WGS84), while you did in the Lambert II projection ? Or did you spot an error in my code (such as when I used the point2grid() function that I am not very familiar with) ? Also, I am curious about how you generated your grid “from scratch” as you said ? Jul 9 at 20:21
  • 1
    I'm not sure exactly how the points2grid algorithm works, but when I tested, I made a new grid with this method, but with the size and extent parameters that I mentioned above, and in Lambert II. Then I checked the distances in Lambert II as well, and all the points were within 3 meters, more or less, of the grid centroids. I would try it with that kind of simple grid, it might give you better results, since from my observations and test, I know it closely matches that specific grid.
    – FSimardGIS
    Jul 9 at 23:42
  • That is indeed much better. With your handmade grid the mean distance between true centroids and estimated centroids ranges from 0 to 5m (in Lambert II, as in WGS84), with an average offset of 2.17m. Therefore I will not use the points2grid() method after all... Thanks a lot for your help ! Jul 10 at 17:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.