CODEWITHSUNDEEP
rspatial
JD Long
JD Long

Reputation: 60756

Creating an equal distance spatial grid in R

I am trying to make an equal dimension square grid over a given area using R. I want my grid to be 1km x 1km square. I see examples like this which illustrate equal lat/long grids:

Creating a regular polygon grid over a spatial extent, rotated by a given angle

but that's not even size. It seems like I should be able to take the st_make_grid function and create this, but I can't grok how to make the grid 1km x 1km.

https://r-spatial.github.io/sf/reference/st_make_grid.html

I'd like, for example, to start at (37,-89.2) and end at (36.2,-86.8) and create an evenly spaced grid that's 1km x 1km. How would I do that with R?

Note: the tricky part it seems, is keeping the grid really 1km x 1km over a very large area. I can keep the grid equal dimensions in decimal degrees, but that's not equal distance on the ground.

I've been able to do this with PostGIS, thanks to a crafty answer here. In PostGIS I've created a function: 

CREATE OR REPLACE FUNCTION public.makegrid_2d (
  bound_polygon public.geometry,
  width_step integer,
  height_step integer
)
RETURNS public.geometry AS
$body$
DECLARE
  Xmin DOUBLE PRECISION;
  Xmax DOUBLE PRECISION;
  Ymax DOUBLE PRECISION;
  X DOUBLE PRECISION;
  Y DOUBLE PRECISION;
  NextX DOUBLE PRECISION;
  NextY DOUBLE PRECISION;
  CPoint public.geometry;
  sectors public.geometry[];
  i INTEGER;
  SRID INTEGER;
BEGIN
  Xmin := ST_XMin(bound_polygon);
  Xmax := ST_XMax(bound_polygon);
  Ymax := ST_YMax(bound_polygon);
  SRID := ST_SRID(bound_polygon);

  Y := ST_YMin(bound_polygon); --current sector's corner coordinate
  i := -1;
  <<yloop>>
  LOOP
    IF (Y > Ymax) THEN  
        EXIT;
    END IF;

    X := Xmin;
    <<xloop>>
    LOOP
      IF (X > Xmax) THEN
          EXIT;
      END IF;

      CPoint := ST_SetSRID(ST_MakePoint(X, Y), SRID);
      NextX := ST_X(ST_Project(CPoint, $2, radians(90))::geometry);
      NextY := ST_Y(ST_Project(CPoint, $3, radians(0))::geometry);

      i := i + 1;
      sectors[i] := ST_MakeEnvelope(X, Y, NextX, NextY, SRID);

      X := NextX;
    END LOOP xloop;
    CPoint := ST_SetSRID(ST_MakePoint(X, Y), SRID);
    NextY := ST_Y(ST_Project(CPoint, $3, radians(0))::geometry);
    Y := NextY;
  END LOOP yloop;

  RETURN ST_Collect(sectors);
END;
$body$
LANGUAGE 'plpgsql';

Then I can call it and pass it a polygon:

SELECT (
    ST_Dump(
      makegrid_2d(
        ST_GeomFromText(
          'Polygon((-75 42, -75 40, -73 40, -73 42, -75 42))',
          4326
        )  ,
         1000, -- width step in meters
         1000  -- height step in meters
       ) 
    )
  ) .geom AS cell into test_grid_cell;

But as you can see, even with PostGIS this is hardly a canned routine. With a bit of effort I think I could port this to sf, but I'm not overly excited about it...

Upvotes: 4

Views: 7938

Answers (1)

Jindra Lacko
Jindra Lacko

Reputation: 8699

Consider this example; it uses the boundaries of the City of Prague for basis of the grid.

The key part is to make sure your sf object is in a metric CRS (or a feety one if you are an American and feeling patriotic); it does not matter which as long as it is in a projected CRS (i.e. st_is_longlat(x) returns FALSE).

library(dplyr)
library(RCzechia) # a package of Czech administrative areas
library(sf)


mesto <- kraje() %>% # All Czech NUTS3 ...
  filter(NAZ_CZNUTS3 == 'Hlavní město Praha') %>% # ... city of Prague
  st_transform(5514) # a metric CRS 

grid_spacing <- 1000  # size of squares, in units of the CRS (i.e. meters for 5514)

polygony <- st_make_grid(mesto, square = T, cellsize = c(grid_spacing, grid_spacing)) %>% # the grid, covering bounding box
  st_sf() # not really required, but makes the grid nicer to work with later

plot(polygony, col = 'white')
plot(st_geometry(mesto), add = T)

gridded Prague

Upvotes: 16

Related Questions