Superposition de polygone spatial avec grille et vérification de la localisation des coordonnées spécifiques à un élément de la grille à l'aide de R [fermé]

32

Comment peut-on utiliser R pour

  1. diviser un fichier de formes en carrés / sous-polygones de 200 mètres,
  2. tracer cette grille (y compris les numéros d’identification pour chaque carré) sur la carte originale ci-dessous, et
  3. Évaluer dans quel carré se trouvent les coordonnées géographiques spécifiques .

Je suis un débutant en système d’information géographique et c’est peut-être une question fondamentale, mais je n’ai pas trouvé de didacticiel expliquant comment procéder dans R.

Jusqu'à présent, j'ai chargé un fichier de formes de NYC et tracé des coordonnées géographiques exemplaires.

Je cherche un exemple (code R) comment procéder avec les données ci-dessous.

# 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")

entrez la description de l'image ici

majom
la source
Voir aussi stackoverflow.com/q/17801398/287948
Peter Krauss Le

Réponses:

36

Voici un exemple utilisant un SpatialGridobjet:

### 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
# ...

Vous pouvez maintenant utiliser la overméthode implémentée pour obtenir les identifiants de cellule:

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

Pour tracer le fichier de formes et la grille avec les ID de cellule:

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 sans couleur / clé de couleur:

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
la source
Cela ressemble à une réponse pour moi, mais au cas où vous cherchiez quelque chose de différent. Essayez la balise r dans stackoverflow stackoverflow.com/search?q=R+tag
Brad Nesom
@rcs ce code ressemble à ce que j'essaye de faire, mais mon fichier de formes a une projection différente: proj4string (DK_reg1) [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" quelqu'un a-t-il des suggestions sur la manière de fractionner ces fichiers de formes de cette projection en 1000 cellules de grille de taille égale? puis sélectionnez au hasard 100 d'entre eux et mettez-les en surbrillance?
Je Del Toro
9

L'ensemble de données de New York fourni dans la question n'est plus disponible au téléchargement. J'utilise l'ensemble de données nc du paquet sf pour illustrer une solution utilisant le paquet 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)

entrez la description de l'image ici

sebdalgarno
la source
Merci. J'ai mis à jour le lien dans ma question pour modifier les modifications sur leur page Web. Maintenant, cela devrait fonctionner à nouveau.
majom
J'ai vraiment besoin de commencer à utiliser le sfpaquet. C'est génial!
philiporlando
Existe-t-il un moyen simple de tracer uniquement les cellules de la grille qui se croisent avec le polygone d'état?
philiporlando
st_intersection (grid_50, nc) devrait le faire
sebdalgarno
Y at-il un moyen de reproduire la même chose, mais les points au centre de chaque grille, donc une grille est dessinée avec la latitude / longitude comme centre de la grille @sebdalgarno
Vijay Ramesh
2

Si vous n'avez pas examiné le package raster R, celui-ci dispose d'outils permettant de convertir en objets SIG vectoriels. Vous devez donc être en mesure de: a) créer un raster (grille) avec des cellules de 200x200m et b) le convertir en un ensemble de un identifiant logique quelconque. À partir de là, je regarderais le paquet sp pour aider à intersecter les points et la grille de polygones. Cette page http://cran.r-project.org/web/packages/sp/vignettes/over.pdf pourrait être un bon début. En parcourant la documentation du paquet sp, vous pourrez peut-être commencer par la classe SpatialGrid et simplement ignorer la partie raster.

Cokrzys
la source
-1

L '"univers SIG" est complexe et comporte de nombreuses normes selon lesquelles vos données doivent être conformes. Tous les "outils SIG" interopèrent selon les normes SIG . Toutes les "données SIG sérieuses" aujourd'hui (2014) sont stockées dans une base de données .

La meilleure façon "d'utiliser R" dans un contexte SIG, avec d'autres outils FOSS , est intégrée à SQL. Les meilleurs outils sont PostgreSQL 9.X (voir PL / R ) et PostGIS .


Vous répondez:

  • Pour importer / exporter des fichiers de forme: utilisez shp2pgsqletpgsql2shp .
  • Pour "diviser un fichier de forme en carrés / sous-polygones de 200 mètres": voir ST_SnapToGrid(),ST_AsRaster() etc. Nous avons besoin de mieux comprendre vos besoins d'exprimer dans une « recette ».
  • vous dites que le besoin "les coordonnées géographiques sont situées" ... peut-être ST_Centroid()des carrés (?) ... Vous pouvez exprimer "plus mathématiquement" pour que je comprenne.

... Peut-être n'avez-vous pas besoin d'une conversion de trame, mais seulement d'une matrice de points échantillonnés régulièrement.


Une méthode primitive consiste à utiliser R sans PL / R , dans un compilateur externe habituel: convertissez uniquement vos polygones et exportez-les sous forme ou sous forme de WKT (voir ST_AsText), puis convertissez les données avec awk ou un autre filtre au format R.

Peter Krauss
la source
1
Merci de votre aide. Cependant, je préférerais fortement une solution qui repose entièrement sur R et les packages existants. Lorsque je suis en mesure de scinder le fichier de forme en point.in.polygonsous-polygones de 200 m * 200 m, je peux vérifier avec quelles coordonnées se trouvent quels polygones. Mon problème est de scinder le fichier de formes original en ces sous-polygones.
Majom