Melapisi poligon spasial dengan kisi dan memeriksa di mana koordinat spesifik elemen kisi berada menggunakan R [ditutup]

Bagaimana seseorang bisa menggunakan R to

1. membagi shapefile dalam kotak / sub-poligon 200 meter,
2. plot kisi ini (termasuk nomor ID untuk setiap kotak) di atas peta asli di bawah, dan
3. mengevaluasi di mana koordinat geografis spesifik persegi berada .

Saya seorang pemula dalam GIS dan ini mungkin pertanyaan mendasar, tapi saya belum menemukan tutorial tentang cara melakukan ini di R.

Apa yang saya lakukan sejauh ini adalah memuat shapefile NYC dan merencanakan beberapa koordinat geografis yang patut dicontoh.

Saya mencari contoh (kode R) bagaimana ini dengan data di bawah ini.

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

Berikut ini contoh menggunakan SpatialGridobjek:

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

Sekarang Anda dapat menggunakan over-method yang diimplementasikan untuk mendapatkan ID sel:

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

Untuk memplot shapefile dan kisi-kisi dengan ID sel:

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

atau tanpa warna / tombol warna:

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

Dataset New York yang disediakan dalam pertanyaan tidak lagi tersedia untuk diunduh. Saya menggunakan dataset nc dari paket sf untuk menunjukkan solusi menggunakan paket 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)

Jika Anda belum melihat paket R raster, ia memiliki alat untuk mengonversi ke / dari objek GIS vektor sehingga Anda harus dapat a) membuat raster (kisi) dengan sel 200x200m dan b) mengonversinya menjadi satu set poligon dengan id logis dari beberapa jenis. Dari sana saya akan melihat paket sp untuk membantu memotong titik dan kisi poligon. Halaman http://cran.r-project.org/web/packages/sp/vignettes/over.pdf ini mungkin merupakan awal yang baik. Berkeliaran melalui dokumen paket sp Anda mungkin dapat mulai dengan kelas SpatialGrid dan hanya melewatkan bagian raster sepenuhnya.

"Alam semesta SIG" rumit dan memiliki banyak standar yang harus sesuai dengan data Anda. Semua "alat GIS" beroperasi dengan standar GIS . Semua "data GIS serius" hari ini (2014) disimpan dalam database .

Cara terbaik untuk "menggunakan R" dalam konteks GIS, dengan alat FOSS lainnya , tertanam ke dalam SQL. Alat terbaik adalah PostgreSQL 9.X (lihat PL / R ) dan PostGIS .

Anda menjawab:

• Untuk mengimpor / mengekspor file bentuk: gunakan shp2pgsqldanpgsql2shp .
• Untuk "membagi file bentuk di 200 meter persegi / sub-poligon": lihat ST_SnapToGrid(), ST_AsRaster(), dll Kami perlu memahami lebih baik kebutuhan Anda untuk mengekspresikan menjadi "resep".
• Anda mengatakan bahwa perlu "koordinat geografis terletak" .. mungkin ST_Centroid()dari kotak (?) ... Anda dapat mengekspresikan "lebih matematis" jadi saya mengerti.

... Mungkin Anda tidak memerlukan konversi raster, hanya matriks poin regurlar-sampel.

Cara primitif adalah menggunakan R tanpa PL / R , dalam kompiler eksternal Anda yang biasa: hanya konversi poligon Anda dan ekspor sebagai bentuk atau sebagai WKT (lihat ST_AsText), lalu konversi data dengan awk atau filter lain ke format R.