Sobreposição de polígono espacial com grade e verificação de quais coordenadas específicas do elemento de grade estão localizadas usando R [fechado]

32

Como se pode usar R para

  1. divida um arquivo de forma em quadrados / sub-polígonos de 200 metros,
  2. plote esta grade (incluindo números de identificação para cada quadrado) sobre o mapa original abaixo e
  3. avaliar em que coordenadas geográficas específicas quadradas estão localizadas .

Sou iniciante em SIG e essa talvez seja uma pergunta básica, mas não encontrei um tutorial sobre como fazer isso em R.

O que fiz até agora foi carregar um arquivo de forma de NYC e plotar algumas coordenadas geográficas exemplares.

Estou procurando um exemplo (código R) de como fazer isso com os dados abaixo.

# Load packages 
library(maptools)

# Download shapefile for NYC
# OLD URL (no longer working)
# shpurl <- "http://www.nyc.gov/html/dcp/download/bytes/nybb_13a.zip"
shpurl <- "https://www1.nyc.gov/assets/planning/download/zip/data-maps/open-data/nybb_13a.zip"

tmp    <- tempfile(fileext=".zip")
download.file(shpurl, destfile=tmp)
files <- unzip(tmp, exdir=getwd())

# Load & plot shapefile
shp <- readShapePoly(files[grep(".shp$", files)])
plot(shp)

# Define coordinates 
points_of_interest <- data.frame(y=c(919500, 959500, 1019500, 1049500, 1029500, 989500), 
                 x =c(130600, 150600, 180600, 198000, 248000, 218000),
                 id  =c("A"), stringsAsFactors=F)

# Plot coordinates
points(points_of_interest$y, points_of_interest$x, pch=19, col="red")

insira a descrição da imagem aqui

majestade
fonte
Veja também stackoverflow.com/q/17801398/287948
Peter Krauss

Respostas:

36

Aqui está um exemplo usando um SpatialGridobjeto:

### read shapefile
library("rgdal")
shp <- readOGR("nybb_13a", "nybb")

proj4string(shp)  # units us-ft
# [1] "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 
# +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83
# +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0"

### define coordinates and convert to SpatialPointsDataFrame
poi <- data.frame(x=c(919500, 959500, 1019500, 1049500, 1029500, 989500),
                  y=c(130600, 150600, 180600, 198000, 248000, 218000),
                  id="A", stringsAsFactors=F)
coordinates(poi) <- ~ x + y
proj4string(poi) <- proj4string(shp)

### define SpatialGrid object
bb <- bbox(shp)
cs <- c(3.28084, 3.28084)*6000  # cell size 6km x 6km (for illustration)
                                # 1 ft = 3.28084 m
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)
grd
# cellcentre.offset 923018 129964
# cellsize           19685  19685
# cells.dim              8      8

sp_grd <- SpatialGridDataFrame(grd,
                               data=data.frame(id=1:prod(cd)),
                               proj4string=CRS(proj4string(shp)))
summary(sp_grd)
# Object of class SpatialGridDataFrame
# Coordinates:
#      min     max
# x 913175 1070655
# y 120122  277602
# Is projected: TRUE
# ...

Agora você pode usar o overmétodo implementado para obter os IDs de célula:

over(poi, sp_grd)
#   id
# 1 57
# 2 51
# 3 38
# 4 39
# 5 14
# 6 28

Para plotar o shapefile e a grade com os IDs de célula:

library("lattice")
spplot(sp_grd, "id",
       panel = function(...) {
         panel.gridplot(..., border="black")
         sp.polygons(shp)
         sp.points(poi, cex=1.5)
         panel.text(...)
       })

spplot1

ou sem cor / chave de cor:

library("lattice")
spplot(sp_grd, "id", colorkey=FALSE,
       panel = function(...) {
         panel.gridplot(..., border="black", col.regions="white")
         sp.polygons(shp)
         sp.points(poi, cex=1.5)
         panel.text(..., col="red")
       })

spplot2

rcs
fonte
Isso me parece uma resposta, mas caso você esteja procurando algo diferente. Experimente o tag r em stackoverflow stackoverflow.com/search?q=R+tag
Brad Nesom
@rcs esse código se parece exatamente com o que estou tentando fazer, mas meu shapefile está em uma projeção diferente: proj4string (DK_reg1) [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" alguém tem alguma sugestão sobre como quebrar esses shapefiles dessa projeção em 1000 células de grade de tamanho igual? e selecione aleatoriamente 100 deles e destaque-os?
I Del Toro
9

O conjunto de dados de Nova York fornecido na pergunta não está mais disponível para download. Eu uso o conjunto de dados nc do pacote sf para demonstrar uma solução usando o pacote sf:

library(sf)
library(ggplot2)

# read nc polygon data and transform to UTM 
nc <- st_read(system.file('shape/nc.shp', package = 'sf')) %>%
  st_transform(32617)

# random sample of 5 points
pts <- st_sample(nc, size = 5) %>% st_sf

# create 50km grid - here you can substitute 200 for 50000
grid_50 <- st_make_grid(nc, cellsize = c(50000, 50000)) %>% 
  st_sf(grid_id = 1:length(.))

# create labels for each grid_id
grid_lab <- st_centroid(grid_50) %>% cbind(st_coordinates(.))

# view the sampled points, polygons and grid
ggplot() +
  geom_sf(data = nc, fill = 'white', lwd = 0.05) +
  geom_sf(data = pts, color = 'red', size = 1.7) + 
  geom_sf(data = grid_50, fill = 'transparent', lwd = 0.3) +
  geom_text(data = grid_lab, aes(x = X, y = Y, label = grid_id), size = 2) +
  coord_sf(datum = NA)  +
  labs(x = "") +
  labs(y = "")

# which grid square is each point in?
pts %>% st_join(grid_50, join = st_intersects) %>% as.data.frame

#>   grid_id                 geometry
#> 1      55 POINT (359040.7 3925435)
#> 2      96   POINT (717024 4007464)
#> 3      91 POINT (478906.6 4037308)
#> 4      40 POINT (449671.6 3901418)
#> 5      30 POINT (808971.4 3830231)

insira a descrição da imagem aqui

sebdalgarno
fonte
Obrigado. Atualizei o link na minha pergunta para refazer as alterações na página da web. Agora deve funcionar novamente.
majom 5/0318
Eu realmente preciso começar a usar o sfpacote. Isso é incrível!
philiporlando
Existe uma maneira fácil de plotar apenas as células da grade que se cruzam com o polígono do estado?
23418 philiporlando
st_intersection (grid_50, nc) deve fazê-lo
sebdalgarno
Existe uma maneira de replicar o mesmo, mas os pontos no centro de cada grade, então uma grade está sendo desenhada com o lat / comprimento como o centro da grade @sebdalgarno
Vijay Ramesh
2

Se você não examinou o pacote raster r, ele possui ferramentas para converter objetos GIS de / para vetor, para que você possa: a) criar uma raster (grade) com células de 200x200 m eb) convertê-la em um conjunto de polígonos com um id lógico de algum tipo. A partir daí, examinaria o pacote sp para ajudar na interseção dos pontos e da grade de polígonos. Esta http://cran.r-project.org/web/packages/sp/vignettes/over.pdf pode ser um bom começo. Percorrendo os documentos do pacote sp, você pode começar com a classe SpatialGrid e pular completamente a parte raster.

cokrzys
fonte
-1

O "universo GIS" é complexo e possui muitos padrões para que seus dados sejam compatíveis. Todas as "ferramentas GIS" interoperam pelos padrões GIS . Todos os "dados GIS sérios" hoje (2014) são armazenados em um banco de dados .

A melhor maneira de "usar R" em um contexto GIS, com outras ferramentas de software livre , é incorporada ao SQL. As melhores ferramentas são o PostgreSQL 9.X (veja PL / R ) e o PostGIS .


Sua resposta:

  • Para importar / exportar arquivos de forma: use shp2pgsqlepgsql2shp .
  • Para "dividir um arquivo de formato em 200 praças metro / sub-polígonos": ver ST_SnapToGrid(), ST_AsRaster()etc. Precisamos entender melhor suas necessidades de expressar em uma "receita".
  • você diz que precisa de "coordenadas geográficas localizadas" ... talvez ST_Centroid()dos quadrados (?) ... Você pode expressar "mais matematicamente" para que eu entenda.

... Talvez você não precise de nenhuma conversão raster, apenas uma matriz de pontos com amostra regular.


Uma maneira primitiva é usar R sem PL / R , em seu compilador externo habitual: converta apenas seus polígonos e exporte como forma ou WKT (consulte ST_AsText), depois converta dados com awk ou outro filtro para o formato R.

Peter Krauss
fonte
1
Obrigado pela ajuda. No entanto, eu preferiria fortemente uma solução que se baseie completamente no R e nos pacotes existentes. Quando consigo dividir o arquivo de formas em subpolígonos de 200m * 200m, posso verificar com point.in.polygonquais coordenadas estão em quais polígonos. Meu problema é dividir o shapefile original nesses sub-polígonos.
majom