От шейп файл с полигони/площи и точки (ширина, дължина), разберете към кой полигон/площ принадлежи всяка точка? В Р

Опитвам се да идентифицирам в кой многоъгълник (ZCTA... известен още като аналог на пощенския код) принадлежи дадена точка, като се има предвид набор от точки и шейпфайл. Въпреки че има няколко въпроса от този тип, почти всички изглежда ме насочват към QGIS. Въпреки че ще отида и ще науча друг инструмент, ако е необходимо, има ли лесен начин да направя това в R? Имам опит в R среда... не толкова в GIS пространството.

Форматът, който използвам, се намира тук: ftp: //ftp.gisdata.mn.gov/pub/gdrs/data/pub/us_mn_state_mngeo/bdry_zip_code_tabulation_areas/shp_bdry_zip_code_tabulation_areas.zip

Първият ми опит беше да заредя шейп файла като SpatialPolygonsDataFrame, точките като SpatialPointsDataFrame, след което да използвам "over()", за да получа индексите на полигоните, които съвпадат:

library(maptools)
library(maps)
library(sp)

mn.zip.map <- readShapePoly("zip_code_tabulation_areas.shp")
# The shapefile is the one referenced in the link above

latlon <- data.frame(matrix(0,nrow=2,ncol=1))
latlon$lat <- c(44.730178, 44.784711)
latlon$lon <- c(-93.235381, -93.476415)
latlon[1] <- NULL
coordinates(latlon) = ~lon+lat
indices <- over(latlon, mn.zip.map)

С резултати:

> indices
ZCTA5CE10 GEOID10 CLASSFP10 MTFCC10 FUNCSTAT10 ALAND10 AWATER10 INTPTLAT10 INTPTLON10
1      <NA>    <NA>      <NA>    <NA>       <NA>      NA       NA       <NA>       <NA>
2      <NA>    <NA>      <NA>    <NA>       <NA>      NA       NA       <NA>       <NA>
      Shape_Leng Shape_Area
1         NA         NA
2         NA         NA

Надявах се да имам изход на първия ред ZCTA5CE10 == 55124 и изход на втория ред ZCTA5CE10 == 55379. Въпреки това, очевидно това не се случва.

Изглежда, че координатните системи не са подравнени... но и двете трябва да са Lat/Lon, нали?

Какво ми липсва тук? Благодаря предварително.


r gis shapefile
person rucker    schedule 13.11.2015    source източник

Отговори (1)


arrow_upward
4
arrow_downward

Мисля, че трябва да настроите и коригирате проекцията:

library(rgdal)
proj4string(mn.zip.map) <- CRS("+proj=utm +zone=15 +datum=NAD83")
mn.zip.map <- spTransform(mn.zip.map, CRS("+proj=longlat"))
proj4string(latlon) <- CRS(proj4string(mn.zip.map))
over(latlon, mn.zip.map)
#   ZCTA5CE10 GEOID10 CLASSFP10 MTFCC10 FUNCSTAT10   ALAND10 AWATER10  INTPTLAT10   INTPTLON10 Shape_Leng Shape_Area
# 1     55124   55124        B5   G6350          S  43572536  1759018 +44.7394617 -093.1938424   27059.59   45295591
# 2     55379   55379        B5   G6350          S 152635134  6181840 +44.7539755 -093.5146083   86609.93  158696544
person lukeA    schedule 14.11.2015
comment
Благодаря, @lukeA! Само за пълнота... има необходимо извикване на библиотека... библиотека(rgdal)... но след това работи брилянтно. - person rucker; 14.11.2015