¿Elegir el valor correcto para proj4string para la lectura de archivos de forma en R?

14

Tengo un archivo shape de polígonos y otro archivo CSV que contiene una lista de puntos como pares (Lat, Lng).

Quiero verificar para cada par (lat, lng) del archivo CSV en qué polígono cae dentro ...

El archivo de forma se proyecta y el archivo de proyecto se lee así:

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]]

Mi plan es el siguiente:

  1. Lea el archivo de forma usando la readShapePolyfunción en el MapToolspaquete R.
  2. Lea las coordenadas de puntos del archivo CSV en un marco de datos y conviértalo a SpatialPointsDataFrame
  3. Use la overfunción para determinar en qué polígono cae dentro.

Para hacerlo, necesito especificar proj4stringmientras se carga el archivo de forma en el paso 1 y también transformar las coordenadas del archivo CSV en el mismo sistema de proyección usando la spTransformfunción antes de aplicar la overfunción en el paso 3, ya que requiere que los puntos y los polígonos deben estar bajo el mismo sistema de proyección.

¿Alguna idea sobre cuál debería ser el valor correcto para el contenido del archivo de proyecto que se muestra arriba?

coordinate-system shapefile r proj Moustafa Alzantot
fuente
Si su archivo de forma ya tiene la proyección definida, use "readOGR" en el paquete rgdal. Este paquete es un contenedor para GDAL y realmente reemplaza la funcionalidad de lectura / escritura de archivos de forma en las herramientas de mapa. Esta función maneja todos los tipos de topología y retiene la información de proyección.
Jeffrey Evans
Cuando intento cargar el archivo de forma usando la readOGRfunción, siempre obtengo el error No se puede abrir el archivo
Moustafa Alzantot
OK, ahora he podido leer el archivo usando readOGR .. el uso de la summaryfunción para el SpatialPolygonDataFrameobjeto me dio el valor correcto para elproj4string
Moustafa Alzantot
Bueno, sin detalles sobre cómo está utilizando la función, ¡es difícil ayudarlo! Parte de la sintaxis es el directorio en el que residen los datos y no necesita la extensión .shp en el nombre del archivo. Algo como readOGR (getwd (), "YourShape") debería funcionar si tiene su directorio de trabajo configurado en el mismo lugar donde está su shepfile.
Jeffrey Evans
Gracias @JeffreyEvans, funcionó ahora y lo usé para obtener el proj4string
Moustafa Alzantot

Respuestas:

14

Proj4string es una cadena PROJ4 crs válida .

ver ¿Cómo puedo obtener la cadena proj4 o el código EPSG de un archivo shapeprile .prj? y Shapefile PRJ a PostGIS SRID tabla de búsqueda?

en breve:

  • Puede usar gdalinfo como en la primera referencia o los enlaces de Python GDAL como en la segunda referencia.

O

  • vaya a Prj2EPSG (un servicio simple para convertir información de proyección de texto conocida de archivos .prj en códigos EPSG estándar)
  • Ingrese el contenido de su archivo prj

ingrese la descripción de la imagen aquí

  • el resultado es EPSG: 27700, por lo que una primera versión de la cadena PROJ4 es

    " + init = epsg: 27700 "

`O

ingrese la descripción de la imagen aquí

  • haga clic en Proj4 y la cadena PROJ4 completa es:

    " + proj = tmerc + lat_0 = 49 + lon_0 = -2 + k = 0.9996012717 + x_0 = 400000 + y_0 = -100000 + ellps = airy + datum = OSGB36 + unidades = m + no_defs "

gene
fuente
10

Aquí hay un sitio web muy útil para recuperar el código EPSG para una proyección dada. En su caso, la proyección es "EPSG: 27700". Si tiene proyecciones definidas para el archivo de forma, puede asignar la proyección al crear el SpatialPointsDataFrame y luego usar la definición de proyección de su archivo de forma importado. El uso de "readOGR" del paquete rgdal retendrá las definiciones de proyección. Aquí hay un ejemplo de asignación y extracción de cadenas de coordenadas en datos de clase 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
fuente