Reproyeccion
Introducción
Los sistemas de referencia de coordenadas (SRC) pueden ser de dos tipos principales: sistemas de coordenadas geográficas (‘lon/lat’, con unidades en grados de longitud y latitud) y proyectadas (normalmente con unidades de metros desde un datum). Este capítulo demuestra cómo establecer y transformar datos geográficos de un CRS a otro y, además, destaca problemas específicos que pueden surgir debido a ignorar CRSs y que deberías tener en cuenta, especialmente si tus datos se almacenan con coordenadas lon/lat.
En muchos proyectos no hay necesidad de preocuparse por los diferentes CRS, y mucho menos de realizar conversiones entre ellos. Sin embargo, si conoces el CRS de tus datos y las consecuencias para las operaciones geométricas (tratadas en la siguiente sección), los CRSs deberían simplemente funcionar entre bastidores: la gente a menudo necesita aprender de repente sobre los CRSs cuando las cosas van mal. Tener un CRS claramente definido en el que estén todos los datos del proyecto, además de entender cómo y por qué usar diferentes CRSs, puede asegurar que las cosas no vayan mal.
Este capítulo enseña los fundamentos de los CRS, demuestra las consecuencias de utilizar diferentes CRS (incluyendo lo que puede salir mal) y cómo “reproyectar” conjuntos de datos de un sistema de coordenadas a otro. En las siguientes secciones introducimos los CRSs en R, cómo obtener y establecer CRSs asociados con objetos espaciales, cuándo reproyectar y qué CRS utilizar, entre otras cosas.
Sistemas de referencia de coordenadas
La mayoría de las herramientas geográficas modernas que requieren conversiones de CRS, incluidos los paquetes básicos de R-spatial y el software SIG de escritorio como QGIS, interactúan con PROJ, una biblioteca C++ de código abierto que “transforma coordenadas de un sistema de referencia de coordenadas (CRS) a otro”. Los CRS pueden describirse de muchas maneras, entre las que se incluyen las siguientes.
- Enunciados simples pero potencialmente ambiguos como “está en coordenadas lon/lat”.
- “proj4 strings” formalizadas, aunque ya obsoletas, como
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs. - Con una cadena de texto identificativa “authority:code” como
EPSG:4326.
Cada una de ellas hace referencia a lo mismo: el sistema de coordenadas “WGS84” que constituye la base de las coordenadas del Sistema de Posicionamiento Global (GPS) y de muchos otros conjuntos de datos. Pero, ¿cuál es el correcto?
La respuesta corta es que la tercera forma de identificar los CRS es preferible: Los paquetes sf (y por extensión stars) y terra, además de muchos otros proyectos de software para trabajar con datos geográficos, entienden EPSG:4326. Además, es corto, fácil de recordar y muy ‘localizable’ en línea. sf entiende el identificador más conciso 4326, pero recomendamos la representación más explícita AUTHORITY:CODE para evitar ambigüedades y proporcionar contexto.
La respuesta más larga es que ninguna de las tres descripciones es suficiente y se necesitan más detalles para manejar y transformar los CRS sin ambigüedades: debido a la complejidad de los SIR, no es posible capturar toda la información relevante sobre ellos en cadenas de texto tan cortas. Por este motivo, el Open Geospatial Consortium desarrolló un formato estándar abierto que se denomina WKT (Well-Known Text). Esto se detalla en un documento de más de 100 páginas que “define la estructura y el contenido de una implementación de cadena de texto del modelo abstracto para sistemas de referencia de coordenadas descrito en ISO 19111:2019”. La representación WKT del SIR WGS84, que tiene el identificador EPSG:4326 es la siguiente:
st_crs("EPSG:4326")
## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]La convención de referirse a los identificadores CRSs de la forma AUTORIDAD:CÓDIGO, que también es utilizada por el software geográfico escrito en otros lenguajes, permite referirse a una amplia gama de sistemas de coordenadas formalmente definidos.1 La autoridad más utilizada en los identificadores CRS es EPSG, acrónimo de European Petroleum Survey Group, que publicó una lista normalizada de CRS (el EPSG fue absorbido por el organismo de petróleo y gas el Comité de Geomática de la Asociación Internacional de Productores de Petróleo y Gas en 2005). Pueden utilizarse otras autoridades en los identificadores SIR. ESRI:54030, por ejemplo, se refiere a la implementación de ESRI de la proyección Robinson, que tiene la siguiente cadena WKT (sólo se muestran las 8 primeras líneas):
Las cadenas WKT son exhaustivas, detalladas y precisas, y permiten almacenar y transformar los SIR sin ambigüedades. Contienen toda la información relevante sobre cualquier SIR, incluido su datum y elipsoide, meridiano principal, proyección y unidades. Antes de la aparición de las definiciones WKT CRS, proj-string era la forma estándar de especificar operaciones de coordenadas y almacenar CRSs. Estas representaciones de cadena, construidas sobre una forma clave=valor (por ejemplo, +proj=longlat +datum=WGS84 +no_defs), ya han sido, o deberían ser en el futuro, sustituidas por representaciones WKT en la mayoría de los casos.
Las versiones recientes de PROJ (6+) todavía permiten el uso de proj-strings para definir operaciones de coordenadas, pero algunas claves proj-string (+nadgrids, +towgs84, +k, +init=epsg:) ya no se soportan o se desaconsejan.
Consulta y establecimiento de sistemas de coordenadas
Veamos cómo se almacenan los CRS en los objetos espaciales de R y cómo pueden consultarse y establecerse. Primero veremos cómo obtener y establecer CRS en objetos geográficos vector, comenzando con el siguiente ejemplo:
El CRS puede recuperarse con la función sf st_crs().
La salida es una lista que contiene dos componentes principales:
-
User input(en este casoWGS 84, un sinónimo deEPSG:4326que en este caso se tomó del fichero de entrada), correspondiente a los identificadores CRS descritos anteriormente -
wkt, que contiene la cadena WKT completa con toda la información pertinente sobre el CRS.
El elemento input es flexible, y dependiendo del fichero de entrada o de la entrada del usuario, puede contener la representación AUTHORITY:CODE (por ejemplo, EPSG:4326), el nombre del CRS (por ejemplo, WGS 84), o incluso la definición proj-string. El elemento wkt almacena la representación WKT, que se utiliza al guardar el objeto en un archivo o al realizar cualquier operación de coordenadas. Arriba, podemos ver que el objeto new_vector tiene el elipsoide WGS84, utiliza el meridiano de Greenwich, y el orden de los ejes de latitud y longitud. En este caso, también tenemos algunos elementos adicionales, como USAGE que explica el área adecuada para el uso de este CRS, e ID que señala el identificador del CRS: EPSG:4326.
La función st_crs también tiene una característica útil: podemos recuperar información adicional sobre el CRS utilizado. Por ejemplo, intente ejecutar
-
st_crs(new_vector)$IsGeographicpara comprobar si el CRS es geográfico o no -
st_crs(new_vector)$units_gdalpara averiguar las unidades del CRS -
st_crs(new_vector)$sridextrae su identificador “SRID” (si está disponible) -
st_crs(new_vector)$proj4stringextrae la representación proj-string
En los casos en los que falta un sistema de referencia de coordenadas (CRS) o se establece un CRS incorrecto, se puede utilizar la función st_set_crs():
La obtención y configuración de los CRS funciona de forma similar para los objetos ráster. La función crs() del paquete terra accede a la información CRS de un objeto SpatRaster (nótese el uso de la función cat() para imprimirla correctamente):
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
my_rast = rast(raster_filepath)
cat(crs(my_rast)) # get CRS
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]La salida es la representación de cadena WKT de CRS. La misma función, crs(), también se puede utilizar para establecer un CRS para objetos raster.
Es importante destacar que las funciones st_crs() y crs() no alteran los valores de las coordenadas ni las geometrías. Su función es únicamente establecer una información de metadatos sobre el CRS del objeto.
En algunos casos, el CRS de un objeto geográfico es desconocido, como es el caso del conjunto de datos london creado en el fragmento de código siguiente:
La salida NA muestra que sf no sabe cuál es el CRS y no está dispuesto a adivinarlo (NA significa literalmente ‘no disponible’). A menos que un CRS se especifique manualmente o se cargue desde una fuente que tenga metadatos CRS, sf no hace ninguna suposición explícita sobre qué sistemas de coordenadas, aparte de decir “no lo sé”. Este comportamiento tiene sentido dada la diversidad de CRS disponibles, pero difiere de algunos enfoques, como la especificación del formato de archivo GeoJSON, que hace la suposición simplificadora de que todas las coordenadas tienen un CRS lon/lat: EPSG:4326. Los conjuntos de datos sin un CRS especificado pueden causar problemas: todas las coordenadas geográficas tienen un sistema de coordenadas y el software sólo puede tomar decisiones correctas sobre el trazado y las operaciones geométricas si sabe con qué tipo de CRS está trabajando. Por lo tanto, de nuevo, es importante comprobar siempre el CRS de un conjunto de datos y establecerlo si falta.
Operaciones geométricas sobre datos proyectados y no proyectados
Desde la versión 1.0.0 de sf, la capacidad de R para trabajar con conjuntos de datos vectoriales geográficos que tienen CRS lon/lat ha mejorado sustancialmente, gracias a su integración con el motor de geometría esférica S2. GEOS se utiliza siempre para datos proyectados y datos sin CRS; para datos geográficos se utiliza S2 por defecto, pero puede desactivarse con sf::sf_use_s2(FALSE).
Para demostrar la importancia de los SIR, crearemos un búfer de 100 km alrededor del objeto “Londres” de la sección anterior. También crearemos un búfer deliberadamente defectuoso con una “distancia” de 1 grado, que equivale aproximadamente a 100 km (1 grado son unos 111 km en el ecuador).
La primera etapa consiste en crear tres buffers alrededor de los objetos london y london_geo creados anteriormente con distancias límite de 1 grado y 100 km (o 100.000 m, que puede expresarse como 1e5 en notación científica) desde el centro de Londres:
En la primera línea de arriba, sf asume que la entrada es proyectada y genera un resultado que tiene un buffer en unidades de grados, lo cual es problemático, como veremos. En la segunda línea, sf usa silenciosamente el motor de geometría esférica S2 para calcular la extensión del buffer usando el valor por defecto de max_cells = 1000 — fijado a 100 en la línea tres — las consecuencias se harán evidentes en breve (ver ?s2::s2_buffer_cells para más detalles).
Para resaltar el impacto del uso del motor de geometría S2 por sf para sistemas de coordenadas (geográficas) no proyectadas, lo desactivaremos temporalmente con el comando sf_use_s2() (que está activado, TRUE, por defecto), en el trozo de código que sigue. Al igual que london_buff_no_crs, el nuevo objeto london_geo es una abominación geográfica: tiene unidades de grados, lo que no tiene sentido en la gran mayoría de los casos:
El mensaje de advertencia anterior indica problemas al realizar operaciones geométricas planas en datos lon/lat. Cuando las operaciones de geometría esférica están desactivadas, con el comando sf::sf_use_s2(FALSE), los buffers (y otras operaciones geométricas) pueden dar como resultado salidas sin valor porque utilizan unidades de latitud y longitud, un mal sustituto de las unidades adecuadas de distancias como los metros.
La distancia entre dos líneas de longitud, llamadas meridianos, es de unos 111 km en el ecuador (ejecute geosphere::distGeo(c(0, 0), c(1, 0)) para encontrar la distancia exacta). En los polos se reduce a cero. En la latitud de Londres, por ejemplo, los meridianos están a menos de 70 km de distancia.
Las líneas de latitud, por el contrario, son equidistantes entre sí independientemente de la latitud: siempre están separadas unos 111 km, incluso en el ecuador y cerca de los polos.
Para las operaciones que implican distancias, como el buffering, la única forma de garantizar un buen resultado (sin utilizar motores de geometría esférica) es crear una copia proyectada de los datos y ejecutar la operación sobre ella.
El resultado es un nuevo objeto idéntico a london, pero reproyectado sobre un CRS adecuado (la British National Grid, que tiene un código EPSG de 27700 en este caso) que tiene unidades de metros. Podemos verificar que el CRS ha cambiado usando st_crs() como sigue (parte de la salida ha sido reemplazada por ...):
Entre los componentes más destacados de esta descripción se incluyen el código EPSG (EPSG: 27700) y la cadena detallada wkt (de la que sólo se muestran las 5 primeras líneas). El hecho de que las unidades del CRS, descritas en el campo LENGTHUNIT, sean metros (en lugar de grados) nos indica que se trata de un CRS proyectado: st_is_longlat(london_proj) devuelve ahora FALSE y las operaciones geométricas sobre london_proj funcionarán sin advertencia. Las operaciones de buffer sobre london_proj utilizarán GEOS y los resultados se devolverán con las unidades de distancia adecuadas. La siguiente línea de código crea un buffer alrededor de datos proyectados de exactamente 100 km:
Las geometrías de los tres objetos london_buff* que tienen un CRS especificado creados anteriormente (london_buff_s2, london_buff_lonlat y london_buff_projected) creados en los trozos de código anteriores se ilustran en la Figura @ref(fig:crs-buf).
Se ve claramente que los búferes basados en s2 y CRS proyectados correctamente no están “aplastados”, lo que significa que cada parte del límite del búfer es equidistante de Londres. Los resultados generados a partir de CRSs lon/lat cuando no se utiliza s2, ya sea porque la entrada carece de CRS o porque sf_use_s2() está desactivada, están muy distorsionados, con el resultado alargado en el eje norte-sur, lo que pone de manifiesto los peligros de utilizar algoritmos que asumen datos proyectados en entradas lon/lat (como hace GEOS). Sin embargo, los resultados generados con S2 también están distorsionados, aunque de forma menos dramática. Los dos límites de los búferes de la figura izquierda están dentados, aunque esto sólo puede ser evidente o relevante en el caso del límite grueso que representa un búfer creado con el argumento s2 max_cells fijado en 100. La lección es que los resultados obtenidos a partir de datos lon/lat con S2 también están distorsionados, aunque de forma menos dramática. La lección es que los resultados obtenidos a partir de datos lon/lat mediante S2 serán diferentes de los resultados obtenidos utilizando datos proyectados. La diferencia entre los buffers derivados de S2 y los buffers derivados de GEOS sobre datos proyectados se reducen a medida que aumenta el valor de max_cells: el valor “correcto” para este argumento puede depender de muchos factores y el valor por defecto 1000 es un valor por defecto razonable. Al elegir los valores de max_cells, la velocidad de cálculo debe equilibrarse con la resolución de los resultados. En situaciones en las que los límites curvos son ventajosos, puede ser apropiado transformar a un CRS proyectado antes de almacenar (o realizar otras operaciones geométricas).
La importancia de los CRS (principalmente si son proyectados o geográficos) y los impactos de la configuración por defecto de sf de utilizar S2 para los buffers en los datos lon/lat quedan claros en el ejemplo anterior. En las secciones siguientes se profundiza en qué CRS utilizar cuando se necesitan CRS proyectados y en los detalles de la reproyección de objetos vectoriales y ráster.
¿Cuándo reproyectar?
La sección anterior mostraba cómo establecer el CRS manualmente, con st_set_crs(london, "EPSG:4326"). Sin embargo, en las aplicaciones del mundo real, los CRS suelen establecerse automáticamente cuando se leen los datos. En muchos proyectos, la principal tarea relacionada con los CRS es transformar objetos, de un CRS a otro. Pero, ¿cuándo deben transformarse los datos? ¿Y en qué CRS? No hay respuestas claras a estas preguntas y la selección del CRS siempre implica compromisos. Sin embargo, se ofrecen algunos principios generales que pueden ayudarle a decidir.
En primer lugar, merece la pena considerar cuándo transformar. En algunos casos la transformación a un CRS geográfico es esencial, como cuando se publican datos en línea con el paquete leaflet. Otro caso es cuando hay que comparar o combinar dos objetos con diferentes CRS, como cuando intentamos encontrar la distancia entre dos objetos con diferentes CRS:
Para que los objetos london y london_proj sean geográficamente comparables, uno de ellos debe transformarse en el CRS del otro. Pero, ¿qué CRS utilizar? La respuesta depende del contexto: muchos proyectos, especialmente los de cartografía web, requieren salidas en EPSG:4326, en cuyo caso merece la pena transformar el objeto proyectado. Si, por el contrario, el proyecto requiere operaciones de geometría plana en lugar del motor de operaciones de geometría esférica (por ejemplo, para crear buffers con bordes suaves), puede merecer la pena transformar los datos con un CRS geográfico en un objeto equivalente con un CRS proyectado, como la British National Grid (EPSG:27700).
¿Qué CRS utilizar?
La cuestión de qué CRS es delicada, y rara vez hay una respuesta “correcta”: “No existen proyecciones que sirvan para todo, todas implican distorsión cuando se alejan del centro del marco especificado”.
A la hora de seleccionar CRS geográficos, la respuesta suele ser WGS84. Se utiliza no sólo para la cartografía web, sino también porque los conjuntos de datos GPS y miles de conjuntos de datos ráster y vectoriales se proporcionan por defecto en este CRS. WGS84 es el CRS más común en el mundo, por lo que conviene conocer su código EPSG: 4326. Este “número mágico” puede utilizarse para convertir objetos con CRS proyectados inusuales en algo que se entienda ampliamente.
¿Qué ocurre cuando se necesita un CRS proyectado? En algunos casos, no es algo que podamos decidir libremente: “a menudo, la elección de la proyección corre a cargo de un organismo cartográfico público”. Esto significa que, al trabajar con fuentes de datos locales, probablemente sea preferible trabajar con el CRS en el que se proporcionaron los datos, para garantizar la compatibilidad, aunque el CRS oficial no sea el más preciso.
Una proyección utilizada habitualmente por defecto es la Universal Transverse Mercator (UTM), un conjunto de CRS que divide la Tierra en 60 cuñas longitudinales y 20 segmentos latitudinales. La proyección transversal de Mercator utilizada por los CRSs UTM es conforme pero distorsiona áreas y distancias con severidad creciente con la distancia desde el centro de la zona UTM. Por lo tanto, ¡recomendamos utilizar UTM sólo cuando su objetivo sea preservar los ángulos de un área relativamente pequeña!
Casi todos los lugares de la Tierra tienen un código UTM, como “60H”, que se refiere al norte de Nueva Zelanda, donde se inventó R. Los códigos EPSG de UTM van secuencialmente de 32601 a 32660 para los lugares del hemisferio norte y de 32701 a 32760 para los lugares del hemisferio sur.
Para mostrar cómo funciona el sistema, vamos a crear una función, lonlat2UTM() para calcular el código EPSG asociado a cualquier punto del planeta como sigue:
El siguiente comando utiliza esta función para identificar la zona UTM y el código EPSG asociado para Auckland y Londres:
Actualmente, también disponemos de herramientas que nos ayudan a seleccionar un CRS adecuado, entre las que se incluye el paquete crssuggest. La función principal de este paquete, suggest_crs(), toma un objeto espacial con CRS geográfico y devuelve una lista de posibles CRS proyectados que podrían utilizarse para el área dada.2. Otra herramienta útil es una página web https://jjimenezshaw.github.io/crs-explorer/ que enumera los CRS en función de la ubicación y el tipo seleccionados.
En los casos en los que no esté claro cuál es el CRS adecuado, la elección del CRS debe depender de las propiedades que sea más importante preservar en los mapas y análisis posteriores. Todos los CRS son de igual área, equidistantes, conformes (con formas que permanecen inalteradas), o alguna combinación de los mismos. Pueden crearse CRS personalizados con parámetros locales para una región de interés y pueden utilizarse múltiples CRS en proyectos cuando ningún CRS se adapte a todas las tareas. Los “cálculos geodésicos” pueden proporcionar un recurso alternativo si no hay CRS adecuados. Independientemente del CRS proyectado que se utilice, los resultados pueden no ser precisos para geometrías que abarquen cientos de kilómetros.
A la hora de decidirse por un SRI personalizado, recomendamos lo siguiente:
- Una proyección acimutal igual al área de Lambert (LAEA) para una proyección local personalizada (establezca la latitud y la longitud del origen en el centro del área de estudio), que es una proyección igual al área en todas las ubicaciones pero distorsiona las formas más allá de miles de kilómetros.
- Proyecciones acimutales equidistantes (AEQD) para una distancia en línea recta específicamente precisa entre un punto y el punto central de la proyección local.
- Proyecciones cónicas conformes de Lambert (LCC) para regiones que abarcan miles de kilómetros, con el cono ajustado para mantener unas propiedades de distancia y área razonables entre las líneas secantes.
- Proyecciones estereográficas (STERE) para regiones polares, pero teniendo cuidado de no confiar en cálculos de área y distancia a miles de kilómetros del centro.
Un posible enfoque para seleccionar automáticamente un CRS proyectado específico para un conjunto de datos local es crear una proyección equidistante azimutal (AEQD) para el punto central del área de estudio. Esto implica crear un CRS personalizado (sin código EPSG) con unidades de metros basado en el punto central de un conjunto de datos. Tenga en cuenta que este método debe utilizarse con precaución: ningún otro conjunto de datos será compatible con el CRS personalizado creado y es posible que los resultados no sean precisos cuando se utilicen en conjuntos de datos extensos que abarquen cientos de kilómetros.
Los principios descritos en esta sección se aplican tanto a los conjuntos de datos vectoriales como a los ráster. Sin embargo, algunas características de la transformación CRS son exclusivas de cada modelo de datos geográficos.
Reproyección de geometrías vectoriales
La reproyección de vectores consiste en transformar las coordenadas de los puntos que forman los vértices de las rectas y los polígonos.
Ahora que se ha creado una versión transformada de london, utilizando la función sf st_transform(), se puede encontrar la distancia entre las dos representaciones de Londres. Puede sorprender que londres y londres2 estén a poco más de 2 km de distancia.
Las funciones para consultar y reproyectar CRS se muestran a continuación con referencia a cycle_hire_osm, un objeto sf de spData que representa “estaciones de carga” donde se pueden alquilar bicicletas en Londres. Los CRS de los objetos sf pueden consultarse y establecerse con la función st_crs(). La salida se imprime como múltiples líneas de texto que contienen información sobre el sistema de coordenadas:
st_crs(cycle_hire_osm)
## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCS["WGS 84",
## DATUM["WGS_1984",
## SPHEROID["WGS 84",6378137,298.257223563,
## AUTHORITY["EPSG","7030"]],
## AUTHORITY["EPSG","6326"]],
## PRIMEM["Greenwich",0,
## AUTHORITY["EPSG","8901"]],
## UNIT["degree",0.0174532925199433,
## AUTHORITY["EPSG","9122"]],
## AUTHORITY["EPSG","4326"]]Los componentes principales del CRS, User input y wkt, se imprimen como una sola entidad, la salida de st_crs() es de hecho una lista con nombre de la clase crs con dos elementos, cadenas de caracteres simples llamadas input y wkt, como se muestra en la salida del siguiente trozo de código:
Se pueden recuperar elementos adicionales con el operador $, incluyendo Name, proj4string y epsg (ver ?st_crs y el tutorial CRS y tranformación en el sitio web de GDAL para más detalles):
La representación WKT almacenada en el elemento $wkt del objeto crs_lnd es la última fuente de verdad. Esto significa que las salidas del trozo de código anterior son consultas de la representación wkt proporcionada por PROJ, en lugar de atributos inherentes del objeto y su CRS.
Tanto los elementos wkt como User Input del CRS se modifican cuando se transforma el CRS del objeto. En el siguiente fragmento de código, creamos una nueva versión de cycle_hire_osm con un CRS proyectado (sólo se muestran las 4 primeras líneas de la salida CRS por brevedad):
El objeto resultante tiene un nuevo CRS con un código EPSG 27700. Pero, ¿cómo averiguar más detalles sobre este código EPSG, o sobre cualquier código?
La impresión de un objeto espacial en la consola devuelve automáticamente su sistema de referencia de coordenadas. Para acceder a él y modificarlo explícitamente, utilice la función st_crs, por ejemplo, st_crs(ciclo_hire_osm).
Reproyección de geometrías raster
Existen diferencias importantes en la reproyección de vectores y rásters: transformar un objeto vectorial implica cambiar las coordenadas de cada vértice, pero esto no se aplica a los datos raster. Los rásters se componen de celdas rectangulares del mismo tamaño (expresadas en unidades cartográficas, como grados o metros), por lo que suele ser impracticable transformar las coordenadas de los píxeles por separado. La reproyección de rásters implica la creación de un nuevo objeto ráster, a menudo con un número de columnas y filas distinto del original. A continuación, hay que volver a estimar los atributos, lo que permite “rellenar” los nuevos píxeles con los valores adecuados. En otras palabras, la reproyección ráster puede considerarse como dos operaciones espaciales separadas: una reproyección vectorial de la extensión ráster a otro CRS y el cálculo de nuevos valores de píxel mediante remuestreo. Por lo tanto, en la mayoría de los casos en que se utilizan datos ráster y vectoriales, es mejor evitar la reproyección de rásteres y reproyectar vectores en su lugar.
La reproyección de los rásters regulares también se conoce como “warping”. Además, existe una segunda operación similar denominada “transformación”. En lugar de remuestrear todos los valores, deja todos los valores intactos pero vuelve a calcular nuevas coordenadas para cada celda del raster, cambiando la geometría de la cuadrícula. Por ejemplo, podría convertir el ráster de entrada (una malla regular) en una malla curvilínea. La operación de transformación puede realizarse en R utilizando el paquete stars.
El proceso de reproyección raster se realiza con project() del paquete terra. project() toma un objeto geográfico (un conjunto de datos raster en este caso) y alguna representación CRS como segundo argumento. Nota al margen: el segundo argumento también puede ser un objeto ráster existente con un CRS diferente.
Veamos dos ejemplos de transformación ráster: el uso de datos categóricos y continuos. Los datos de cobertura del suelo suelen representarse mediante mapas categóricos. El archivo nlcd.tif proporciona información para una pequeña área en Utah, EE.UU. obtenida de National Land Cover Database 2011 en el CRS NAD83 / UTM zona 12N, como se muestra en la salida del fragmento de código a continuación:
En esta región se distinguieron 8 clases de ocupación del suelo (la lista completa de las clases de ocupación del suelo de la NLCD2011 puede consultarse en mrlc.gov):
Al reproyectar rásters categóricos, los valores estimados deben ser los mismos que los del original. Para ello se puede utilizar el método del vecino más cercano (near), que establece el valor de cada nueva celda en el valor de la celda más cercana (centro) del ráster de entrada. Un ejemplo es la reproyección de cat_raster a WGS84, un CRS geográfico muy adecuado para la cartografía web. El primer paso es obtener la definición PROJ de este CRS, lo que puede hacerse, por ejemplo, utilizando la página web http://spatialreference.org. El último paso es reproyectar el ráster con la función project() que, en el caso de datos categóricos, utiliza el método del vecino más cercano (near):
Muchas propiedades del nuevo objeto difieren del anterior, incluido el número de columnas y filas (y, por tanto, el número de celdas), la resolución (transformada de metros a grados) y la extensión.
| CRS | nrow | ncol | ncell | resolution | unique_categories |
|---|---|---|---|---|---|
| NAD83 | 1359 | 1073 | 1458207 | 31.5275 | 8 |
| WGS84 | 1246 | 1244 | 1550024 | 0.0003 | 9 |
La reproyección de rásters numéricos (con valores numéricos o en este caso enteros) sigue un procedimiento casi idéntico. Esto se demuestra a continuación con srtm.tif de the Shuttle Radar Topography Mission (SRTM), que representa la altura en metros sobre el nivel del mar (elevación) con el WGS84 CRS:
con_raster = rast(system.file("raster/srtm.tif", package = "spDataLarge"))
crs(con_raster)
## [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"Reproyectaremos este conjunto de datos en un CRS proyectado, pero no con el método del vecino más próximo, que es apropiado para datos categóricos. En su lugar, utilizaremos el método bilineal, que calcula el valor de la celda de salida basándose en las cuatro celdas más cercanas de la trama original. Los valores del conjunto de datos proyectado son la media ponderada por distancia de los valores de estas cuatro celdas: cuanto más cerca esté la celda de entrada del centro de la celda de salida, mayor será su peso. Los siguientes comandos crean una cadena de texto que representa WGS 84 / UTM zona 12N, y reproyectan el ráster en este CRS, utilizando el método bilineal:
con_raster_ea = project(con_raster, "EPSG:32612", method = "bilinear")
crs(con_raster_ea)
## [1] "PROJCRS[\"WGS 84 / UTM zone 12N\",\n BASEGEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]],\n CONVERSION[\"UTM zone 12N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-111,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n AREA[\"Between 114°W and 108°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Alberta; Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA).\"],\n BBOX[0,-114,84,-108]],\n ID[\"EPSG\",32612]]"La reproyección raster sobre variables numéricas también provoca pequeños cambios en los valores y propiedades espaciales, como el número de celdas, la resolución y la extensión. Estos cambios se muestran en la Tabla @ref(tab:rastercrs)3:
| CRS | nrow | ncol | ncell | resolution | mean |
|---|---|---|---|---|---|
| WGS84 | 457 | 465 | 212505 | 0.0008 | 1842.548 |
| UTM zone 12N | 515 | 422 | 217330 | 83.5334 | 1842.083 |
Por supuesto, las limitaciones de las proyecciones de la Tierra en 2D se aplican tanto a los datos vectoriales como a los ráster. En el mejor de los casos, podemos cumplir dos de las tres propiedades espaciales (distancia, área, dirección). Por lo tanto, la tarea a realizar determina qué proyección elegir. Por ejemplo, si estamos interesados en una densidad (puntos por celda de cuadrícula o habitantes por celda de cuadrícula) deberíamos utilizar una proyección de igual área.
Notas
Se pueden utilizar otras formas de referirse a CRS únicos, con cinco tipos de identificadores (código EPSG, PostGIS SRID, INTERNAL SRID, proj-string y cadenas WKT) aceptados por QGIS y otros tipos de identificadores como una variante más verbosa del identificador
EPSG:4326,urn:ogc:def:crs:EPSG::4326.↩︎Este paquete también permite averiguar el verdadero CRS de los datos sin ninguna información CRS adjunta↩︎
Otro cambio menor es que la clase de los valores en el nuevo conjunto de datos ráster proyectado es
numérico. Esto se debe a que el métodobilinealtrabaja con datos continuos y los resultados raramente se convierten en valores enteros. Esto puede tener implicaciones en el tamaño de los archivos cuando se guardan los conjuntos de datos ráster.↩︎