Создание равноправной пространственной сетки в R - PullRequest
Новая
       18

Создание равноправной пространственной сетки в R

0 голосов
/ 15 декабря 2018

Я пытаюсь создать квадратную сетку равного размера на заданной области, используя R. Я хочу, чтобы моя сетка была квадратной 1 км х 1 км.Я вижу подобные примеры, которые иллюстрируют равные сетки широты и долготы:

Создание правильной многоугольной сетки в пространственном экстенте, повернутой на заданный угол

, но это даже неразмер.Кажется, что я должен быть в состоянии взять функцию st_make_grid и создать ее, но я не могу понять, как сделать сетку 1 км х 1 км.

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

Я хотел бы, например, начать с (37, -89.2) и заканчиваться в (36.2, -86.8) и создать равномерно распределенную сетку размером 1 км х 1 км.Как бы я сделал это с R?

Примечание: кажется, что хитрая часть - это то, что сетка действительно 1 км х 1 км на очень большой площади.Я могу держать сетку равными размерами в десятичных градусах, но это не равное расстояние на земле.

Я смог сделать это с PostGIS, благодаря хитрому ответу . В PostGIS я создал функцию: 

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';

Затем яможно назвать его и передать многоугольник: 

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;

Но, как вы можете видеть, даже с PostGIS это вряд ли обычная процедура.Приложив немного усилий, я думаю, что смогу перенести это на sf, но я не слишком рад этому ...

1 Ответ

0 голосов
/ 16 декабря 2018

Рассмотрим этот пример;он использует границы города Праги для основы сетки.

Ключевые части:
- убедитесь, что ваш sf объект находится в метрической системе координат (или в случае, еслиамериканец и патриотизм)
- подсчитать необходимое количество квадратов
- запустить код и наслаждаться ... 

library(tidyverse)
library(RCzechia)
library(sf)


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

krabicka <- st_bbox(mesto)

xrange <- krabicka$xmax - krabicka$xmin
yrange <- krabicka$ymax - krabicka$ymin

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

rozmery <- c(xrange/grid_spacing , yrange/grid_spacing) %>% # number of polygons necessary
  ceiling() # rounded up to nearest integer

polygony <- st_make_grid(mesto, square = T, n = rozmery) %>% # the grid, covering bounding box
  st_intersection(mesto) %>% # only the inside part 
  st_sf() # not really required, but makes the grid nicer to work with later

plot(polygony, col = 'white')

gridded prague

...