Wybierając prawidłową wartość dla proj4string do odczytu pliku kształtu w R?

14

Mam plik kształtu wielokątów i inny plik CSV, który zawiera listę punktów w postaci par (Lat, Lng) ..

Chcę sprawdzić dla każdej pary (łac., Lng) z pliku CSV, w którym wielokącie wpada.

Plik shapefile jest wyświetlany, a plik proj wygląda następująco:

PROJCS["Transverse_Mercator",GEOGCS["GCS_OSGB 1936",
DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["Meter",1]]

Mój plan jest następujący:

  1. Odczytaj plik kształtu za pomocą readShapePolyfunkcji w MapToolspakiecie R.
  2. Odczytaj współrzędne punktów z pliku CSV do ramki danych i przekonwertuj ją na SpatialPointsDataFrame
  3. Użyj overfunkcji, aby określić, w którym wielokącie wchodzi.

Aby to zrobić, muszę określić proj4stringpodczas ładowania pliku kształtu w kroku 1, a także przekształcić współrzędne z pliku CSV do tego samego systemu projekcji za pomocą spTransformfunkcji przed zastosowaniem overfunkcji w kroku 3, ponieważ wymaga to, aby punkty i wielokąty musiały być pod tym samym systemem projekcji.

Masz pojęcie o tym, jaka powinna być poprawna wartość dla zawartości pliku proj pokazanej powyżej?

Moustafa Alzantot
źródło
Jeśli twoje pliki kształtu mają już zdefiniowane odwzorowanie, użyj „readOGR” w pakiecie rgdal. Ten pakiet jest opakowaniem dla GDAL i naprawdę zastępuje funkcję odczytu / zapisu pliku kształtu w maptools. Ta funkcja obsługuje wszystkie typy topologii i zachowuje informacje o rzutowaniu.
Jeffrey Evans
Kiedy próbuję loadign plik kształtu za pomocą readOGRfunkcji I a zawsze dostać nie można otworzyć pliku błędu
Moustafa Alzantot
OK, teraz byłem w stanie odczytać plik przy użyciu readOGR .. użycie summaryfunkcji dla SpatialPolygonDataFrameobiektu dało mi poprawną wartość dlaproj4string
Moustafa Alzantot
Cóż, bez szczegółowych informacji na temat korzystania z tej funkcji trudno jest Ci pomóc! Częścią składni jest katalog, w którym znajdują się dane i nie potrzebujesz rozszerzenia .shp w nazwie pliku. Coś takiego jak readOGR (getwd (), „YourShape”) powinno działać, jeśli katalog roboczy jest ustawiony w tym samym miejscu, w którym znajduje się plik shepfile.
Jeffrey Evans
Dzięki @JeffreyEvans, to działało teraz i użyłem go, aby uzyskać proj4string
Moustafa Alzantot

Odpowiedzi:

14

Łańcuch proj4string jest poprawnym ciągiem crs PROJ4 .

zobacz Jak mogę pobrać ciąg proj4 lub kod EPSG z pliku .prj shapefile? i Shapefile PRJ do tabeli odnośników PostGIS SRID?

w skrócie:

  • Możesz użyć gdalinfo jak w pierwszym odwołaniu lub powiązań GDAL Python jak w drugim odwołaniu.

Lub

  • przejdź do Prj2EPSG (prosta usługa do konwersji znanych informacji o wyświetlaniu tekstu z plików .prj na standardowe kody EPSG)
  • Wpisz treść pliku prj

wprowadź opis zdjęcia tutaj

  • wynikiem jest EPSG: 27700, więc pierwsza wersja ciągu PROJ4 to

    + init = epsg: 27700

„Lub

wprowadź opis zdjęcia tutaj

  • kliknij Proj4, a pełny ciąg PROJ4 to:

    + proj = tmerc + lat_0 = 49 + lon_0 = -2 + k = 0,9996012717 + x_0 = 400000 + y_0 = -100000 + ellps = przewiewny + dane = OSGB36 + jednostki = m + no_defs

gen
źródło
10

Oto bardzo przydatna strona internetowa do pobierania kodu EPSG dla danej projekcji. W twoim przypadku projekcja to „EPSG: 27700”. Jeśli masz zdefiniowane rzuty dla pliku kształtu, możesz przypisać rzut podczas tworzenia SpatialPointsDataFrame, a następnie użyć definicji rzutu z importowanego pliku kształtu. Użycie „readOGR” z pakietu rgdal zachowa definicje projekcji. Oto przykład przypisywania i ciągnięcia ciągów współrzędnych na danych klasy sp.

require(sp)
require(rgdal)

# Use meuse dataset
data(meuse)

# Coerce into SpatialPointsDataframe
coordinates(meuse) <- ~x+y

# Assign projection
proj4string(meuse) <- CRS("+init=epsg:28992")

# Pull some observations and transform to Lat/Long
meuse.geo <- meuse[sample(dim(meuse)[1],10),]
  prj.LatLong <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
    meuse.geo <- spTransform(meuse.geo, prj.LatLong)

# Pull projection string from meuse.geo and use in spTransform
#   to reproject meuse to lat/long  
( prj <- proj4string(meuse.geo) )   
meuse <- spTransform(meuse, CRS(prj))   
Jeffrey Evans
źródło