PREPROCESAMIENTO DE IMÁGENES DE SATÉLITE: LAS PROYECCIONES

G174: FOTOINTERPRETACIÓN Y TELEDETECCIÓN

1 INTRODUCCIÓN

En cualquier aspecto del análisis geoespacial, verificar y garantizar el ajuste espacial y la conformidad de las proyecciones es una tarea fundamental. Para ello, debe definirse el sistema de referencia de coordenadas (en inglés CRS).

2 LOS SISTEMAS DE REFERENCIA DE COORDENADAS

Un sistema de referencia de coordenadas (CRS) puede ser tan simple y arbitrario como la cuadrícula creada para muestrear una parcela o tan complejo como para utilizarse sobre todo el planeta. Se puede definir como un marco de referencia local, regional o global, basado en coordenadas, que se utiliza para representar gráficamente (en un mapa) las ubicaciones reales de entidades geográficas.

Los criterios para diferenciar los CRS son diversos, pero esencialmente se puede recurrir al :

• El marco de medición, que puede ser geográfico (coordenadas esféricas) o planimétrico (coordenadas proyectadas sobre un plano).

• La unidad de medición, generalmente metros para sistemas de coordenadas proyectadas o grados (latitud/longitud) para coordenadas geográficas.

Además, los CRS presentan otras características adicionales, como un esferoide de referencia, el datum y otros parámetros de proyección (uno o más paralelos estándar, meridiano central o posibles cambios en las direcciones de x e y) etc…

A partir de esos criterios, los sistemas de referencia de coordenadas se clasifican en dos grandes grupos:

• Sistema de coordenadas geográficas: las ubicaciones se definen en términos de valores de latitud y longitud (coordenadas medidas desde el centro de la Tierra) sobre una superficie esférica, de tres dimensiones, siendo la unidad de medida los grados. Esto supone que distancias y otros parámetros topográficos no pueden evaluarse con precisión.

• Sistemas de coordenadas proyectadas, también denominadas proyecciones cartográficas. Se basan en la proyección de una superficie curva sobre una superficie plana bidimensional, empleando fórmulas matemáticas para relacionar las coordenadas esféricaso con las coordenadas cartesianas (x, y) que definen las ubicaciones de las entidades geográficas. A diferencia del sistema de coordenadas geográficas, el sistema de coordenadas proyectadas permite cálculos como distancias, áreas constantes etc, pero todas incorporan algunas distorsiones en esas magnitudes. Algunas de las proyecciones más famosas son la ‘Mercator’, ‘UTM’, ‘Robinson’, ‘Lambert’, ‘Sinusoidal’ ‘Robinson’ y ‘Albers’. No existe una proyección mejor que otras; algunas pueden usarse para todo el planeta, mientras que otras son apropiadas solo para áreas pequeñas. Una de las características más importantes de la proyección es si es corresponde a las conocidas como de “áreas iguales” (la escala del mapa es constante) o las “conformes” (las formas como se ven en un globo). Ninguna proyección puede ser a la vez conforme y de áreas iguales (pero pueden ser aproximadamente ambas para áreas más pequeñas, por ejemplo, UTM o Lambert Equal Area para un área grande).

Tradicionalmente, las proyecciones en R usaron la notación [PROJ.4[(ftp://ftp.remotesensing.org/proj/OF90-284.pdf ), nombre de una biblioteca de software de código abierto. En esta notación, cada proyección era identificada mediante una cadena de caracteres que almacena los parámetros que describen tanto la proyección como el sistema de coordenadas. Este formato es conocido como proj4string. Por ejemplo, un sistema de coordenadas geográfico (no proyectado) que utiliza la proyección UTM para el centro de la Península Ibérica (WGS 84/UTM 30N) se define por la siguiente cadena de caracteres “+proj=longlat + datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0”.

• proj=longlat: los datos utilizan un sistema de coordenadas geográficas (latitud y longitud).

• datum=WGS84: el datum WGS84 (el datum hace referencia a un punto de origen de los ejes de coordenadas).

• ellps=WGS84: el elipsoide (como se calcula la esfericidad de la Tierra).

Hoy en día, a los sistemas de proyección más utilizados se les asigna un código EPSG, un código corto compilado originalmente por el European Petroleum Survey Group. Los códigos EPSG con los Sistemas de Referencia de Coordenadas (CRS) utilizados en la representación cartográfica de los servicios Web de Mapas (WMS) del Nodo IDE del MAPAMA son los siguientes (pulsar el enlance para ampliar información sobre cada código):

En coordenadas geográficas: EPSG:4230 (ED50) EPSG:4326 (WGS84) EPSG:4258 (ETRS89)

En coordenadas UTM: EPSG:32628 (WGS84 / UTM zone 28N) EPSG:32629 (WGS84 / UTM zone 29N) EPSG:32630 (WGS84 / UTM zone 30N) EPSG:32631 (WGS84 / UTM zone 31N) EPSG:25828 (ETRS89 / UTM zone 28N) EPSG:25829 (ETRS89 / UTM zone 29N) EPSG:25830 (ETRS89 / UTM zone 30N) EPSG:25831 (ETRS89 / UTM zone 31N) EPSG:23028 (ED50 / UTM zone 28N) EPSG:23029 (ED50 / UTM zone 29N) EPSG:23030 (ED50 / UTM zone 30N) EPSG:23031 (ED50 / UTM zone 31N)

Podemos conocer algo más sobre los sistemas de referencia en las siguientes páginas web: Sistemas de referencia. Coordinate Systems Worldwide En qué zona UTM me encuentro

Página web Spatial Reference

Es muy habitual que, a la hora de trabajar que diferentes fuentes de información, los objetos espaciales no estén proyectados, o puedan estarlo en diferentes CRS. Por ello, deben transformarse en un CRS común para realizar la mayoría de las operaciones típicas de un entorno SIG. Los casos más habituales de transformación de CRS son los siguientes:

• La asignación de un nuevo CRS a un objeto que no disponía de CRS previo.

• La transformación del CRS de un objeto a otro sistema de proyección.

3 CRS DE OBJETOS RÁSTER

En esta sesión se trabajará con los paquetes sf y raster.

library(sf)

## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE

library(raster)

## Loading required package: sp

El paquete raster contiene dos funciones relacionadas con las proyecciones y los CRS:

• Para conocer el CRS:

  • projection() • crs()

• Para asignar un CRS

  • projection(x) <- CRS(“+init=epsg:28992”)

• Para cambiar el CRS de un objeto raster en otro CRS:

  • nuevoRaster <- projectRaster(antiguoRaster, proj4string(nuevoCRS))

3.1 Asignación del CRS a un objeto ráster sin proyección.

El ejercicio comienza creado un objeto ráster a partir de una matriz de datos.

matriz <- matrix(1:36, 6, 6)
matriz

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    7   13   19   25   31
## [2,]    2    8   14   20   26   32
## [3,]    3    9   15   21   27   33
## [4,]    4   10   16   22   28   34
## [5,]    5   11   17   23   29   35
## [6,]    6   12   18   24   30   36

borrame <- raster(matriz)

Una vez creada, el objeto matriz puede ser eliminado ya que no será utilizado

rm(matriz)

Se puede consultar el CRS de un objeto ráster con las siguientes funciones.

projection(borrame)

## [1] NA

crs(borrame) 

## Coordinate Reference System:
## Deprecated Proj.4 representation: NA

La información que proporcionan no es exactamente la misma, pero ambas funciones informan que el objeto ráster no posee sistema de referencia alguno. Para asignarle un CRS se puede actuar de varias maneras.

• Una posible sintaxis asigna una determinada proyección al raster utilizando el formato PROJ.4. Como hemos señalado, este se utiliza cada vez menos.

crs(borrame) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" 
crs(borrame)

## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

• Otra posibilidad es usar el nuevo formato EPSG,

  • Opción 1: con la notación "+init=epsg:codigo_EPSG"

projection(borrame) <- CRS("+init=epsg:23030")

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum European_Datum_1950 in Proj4 definition

crs(borrame)

## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=30 +ellps=intl +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["ED50 / UTM zone 30N",
##     BASEGEOGCRS["ED50",
##         DATUM["European Datum 1950",
##             ELLIPSOID["International 1924",6378388,297,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4230]],
##     CONVERSION["UTM zone 30N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]],
##         ID["EPSG",16030]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Europe - between 6°W and 0°W - Channel Islands (Jersey, Guernsey); France offshore; Gibraltar; Ireland offshore; Norway including Svalbard - offshore; Spain - onshore; United Kingdom - UKCS offshore."],
##         BBOX[35.26,-6,80.53,0]]]

• Opción 2: escribiendo directamente el código EPSG.

projection(borrame) <- CRS("EPSG:25830")

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum European Terrestrial Reference System 1989 in
## Proj4 definition

crs(borrame)

## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=30 +ellps=GRS80 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["ETRS89 / UTM zone 30N",
##     BASEGEOGCRS["ETRS89",
##         DATUM["European Terrestrial Reference System 1989",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4258]],
##     CONVERSION["UTM zone 30N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Europe between 6°W and 0°W: Faroe Islands offshore; Ireland - offshore; Jan Mayen - offshore; Norway including Svalbard - offshore; Spain - onshore and offshore."],
##         BBOX[35.26,-6,80.53,0]],
##     ID["EPSG",25830]]

• En la misma línea, también se podría utilizar la sintaxis siguente

crs(borrame) <- CRS("EPSG:23030")

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum European Datum 1950 in Proj4 definition

Cuando se trabaja con un gran número de objetos al mismo tiempo, puede ser conveniente crear un objeto con la proyección a la que se desea modificar todos ellos, y luego asignárselos a cada objeto mediante la función raster::projectRaster. Podrían utilizarse distintas sintaxis:

• Asignando el nuevo CRS mediante la notación proj4string.

proyeccion1 <- "+proj=longlat +datum=WGS84 +no_defs"
borrame_proy1 <- projectRaster(borrame, 
                         crs = proyeccion1)
crs(borrame_proy1)

## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

• Es más cómodo la variante con la notación EPSG.

proyeccion2 <- "EPSG:4326"
borrame_proy2 <- projectRaster(borrame, 
                         crs = proyeccion2)
crs(borrame_proy2)

## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

El paquete raster dispone de la función raster::compareCRS para determinar si las proyecciones de dos objetos espaciales son idénticas.

compareCRS(borrame_proy1, borrame_proy2)

## [1] TRUE

Para no llenar el Global Enviroment de objetos sin uso, se procede a su eliminación.

rm(matriz, borrame, borrame_proy1, borrame_proy2, proyeccion1, proyeccion2)

## Warning in rm(matriz, borrame, borrame_proy1, borrame_proy2, proyeccion1, :
## objeto 'matriz' no encontrado

3.2 Transformación del CRS de un ráster a otro CRS (reproyección)

La proyección y reproyección de datos ráster supone un problema que no tienen los objetos vectoriales: las nuevas ubicaciones de los píxeles no tienen por qué coincidir con las localizaciones antiguas una vez reproyectados, ya que han sido “obligadas” a encajarse en otras cuadrículas con una disposición diferente. Los píxeles de un objeto ráster poseen el mismo tamaño (en términos de las unidades del CRS aunque su tamaño real puede variar), por lo que no es posible transformar celda a celda.

El cálculo de los valores de los píxeles reproyectados se puede realizar a partir de los valores de los píxeles del ráster original. Por lo tanto, será necesario interpolar esos valores. Esto significa calcular un valor promedio ponderado de los píxeles circundantes más cercanos (cuatro píxeles en la práctica). Alternativamente, se podría usar el método de remuestreo del vecino más cercano, que recoge el valor del píxel más cercano, evitando así el promedio, pero esto tiende a dar resultados bastante pixelados y es formalmente menos correcto. Se recomienda:

• Usar como método el del vecino más cercano (argumento “ngb”) si los valores son categóricos (por ejemplo, usos del suelo).

• Si la variable del ráster es continúa (por ejemplo, la altitud) la interpolación se realiza con el método bilineal (argumento “bilinear”).

Además del suavizado de los valores originales, la reproyección de rásteres con cualquier método que no sea el del vecino más cercano es irreversible: no hay forma de volver a los valores originales. Por lo tanto, si tiene que elegir entre volver a proyectar datos vectoriales para que coincidan con su proyección ráster o volver a proyectar el ráster para que coincida con la proyección vectorial, se recomienda optar por proyectar los objetos vectoriales, que es, además, mucho más rápida.

A continuación se diseñará un objeto ráster formado por 81 celdas, a las que se incorporarán las coordenadas geográficas.

raster <- raster(ncol= 15,                # Número de columnas del ráster
               nrow= 15,                  # Número de filas
               xmx= -2,xmn= -5,           # Coordenadas del eje X (considerémoslas la longitud)
               ymx= 42,ymn= 39)           # Coordenadas del eje y (idem latitud)

Este objeto ráster se rellenará con los siguientes valores

values(raster) <- as.integer(runif(ncell(raster),
                                 min = 1,
                                 max= 6))

Llamando al objeto ráster se puede conocer las características de dicho objeto.

raster

## class      : RasterLayer 
## dimensions : 15, 15, 225  (nrow, ncol, ncell)
## resolution : 0.2, 0.2  (x, y)
## extent     : -5, -2, 39, 42  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 5  (min, max)

Por último, se representará gráficamente este objeto.

plot(raster)
text(raster, digits=1)

Para definir el CRS destino se debería usar a partir de este momento el código EPSG. Esto puede realizarse con varios tipos de sintaxis.

• Creación de un objeto R que contiene el nuevo CRS. Esto se hace con el texto con formato +init=epsg:CODIGO_EPSG; la función agreta automáticamente todos los parámetros de dicho CRS (ED50 / UTM zone 30N).

utm <- crs("+init=epsg:23030")                
utm

## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=30 +ellps=intl +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["ED50 / UTM zone 30N",
##     BASEGEOGCRS["ED50",
##         DATUM["European Datum 1950",
##             ELLIPSOID["International 1924",6378388,297,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4230]],
##     CONVERSION["UTM zone 30N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]],
##         ID["EPSG",16030]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Europe - between 6°W and 0°W - Channel Islands (Jersey, Guernsey); France offshore; Gibraltar; Ireland offshore; Norway including Svalbard - offshore; Spain - onshore; United Kingdom - UKCS offshore."],
##         BBOX[35.26,-6,80.53,0]]]

La reproyección se realiza con la función raster::projectRaster, utilizando para ello el objeto utm con el CRS deseado. La reproyección puede realizarse mediante dos métodos: el argumento method = ngb significa que se utiliza el método del vecino más próximo.

raster_proj_ngb <- projectRaster(from = raster, 
                             crs = utm,
                             method= "ngb") 

De la misma manera, el argumento method = bilinear significa que la interpolación sigue el método bilineal.

raster_proj_bil <- projectRaster(from = raster, 
                             crs = utm,
                             method= "bilinear") 

En la siguiente imagen se representa el objeto ráster original y los resultantes de la reproyección con diferentes métodos.

# Dos filas, dos columnas
par(mfcol = c(2, 2))

# Gráficos
plot(raster, main = "Datos originales")
text(raster, digits=0, cex = 0.5)
plot(raster_proj_ngb, main = "Reproyección ngb")
text(raster_proj_ngb, digits=0, cex = 0.5)
plot(raster_proj_bil, main = "Reproyección bil")
text(raster_proj_bil, digits=0, cex = 0.5)

# Volvemos al estado original
par(mfcol = c(1, 1))

Obsérvese que, mientras en los datos originales el CRS sigue coordenadas geográficas (valores de latitud y longitud), en los objetos reproyectados con el nuevo CRS las unidades son metros.

Otra posibilidad de reproyección consiste en crear un objeto ráster vacío y rellenarlo con los datos reproyectados de un fichero original. Esto permite que la alineación del nuevo objeto sea más correcta y precisa. El objeto raster_nuevo se crea a partir del objeto raster original raster, pero como objeto vacío, con la función projectExtent. Esta última función está concebida para crear objetos ráster vacíos, que sirvan para posteriores análisis. Requiere dos argumentos:

• El nombre del objeto ráster cuyas dimensiones se desea “clonar”.

• El CRS que se desea asignar

raster_nuevo <- projectExtent(raster, 
                              utm)
raster_nuevo

## class      : RasterLayer 
## dimensions : 15, 15, 225  (nrow, ncol, ncell)
## resolution : 17302.64, 22328.6  (x, y)
## extent     : 326919.2, 586458.8, 4316989, 4651918  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +ellps=intl +units=m +no_defs

A este objeto vacío se le puede modificar la resolución original,

res(raster_nuevo)

## [1] 17302.64 22328.60

La nueva resolución espacial pasará a ser de 10000 metros en las dos dimensiones x e y; si quisiéramos píxeles con diferentes resolución espacial (x,y) la sintaxis podría ser c(10000,15000).

res(raster_nuevo) <- 10000 
raster_nuevo

## class      : RasterLayer 
## dimensions : 33, 26, 858  (nrow, ncol, ncell)
## resolution : 10000, 10000  (x, y)
## extent     : 326919.2, 586919.2, 4321918, 4651918  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +ellps=intl +units=m +no_defs

A continuación se reproyectará el objeto ráster original con la resolución del raster vacío; sus valores se usarán para rellenar un nuevo objeto ráster raster_nuevo_resol1000

raster_nuevo_resol1000 <- projectRaster(raster, 
                                      raster_nuevo, 
                                      method = 'ngb')
raster_nuevo_resol1000

## class      : RasterLayer 
## dimensions : 33, 26, 858  (nrow, ncol, ncell)
## resolution : 10000, 10000  (x, y)
## extent     : 326919.2, 586919.2, 4321918, 4651918  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +ellps=intl +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 5  (min, max)

Por último, para comprobar el resultado se procederá a su representación gráfica.

plot(raster_nuevo_resol1000, main = "Reproyección")
text(raster_nuevo_resol1000, digits=0, cex = 0.4)

rm(raster_nuevo, raster_nuevo_resol1000, raster_proj_bil, raster_proj_ngb, utm)

4 PROYECCIÓN DE OBJETOS VECTORIALES

Es muy habitual también combinar objetos en formato ráster con objetos en formato vectorial. Es necesario, por lo tanto, que ambos formatos compartan el mismo CRS.

4.1 Creación de un objeto vectorial

Podemos crear un objeto vectorial a partir de un objeto ráster. La función extent() extrae de un rasterLayer o un objeto espacial (vector) un objeto extent, en realidad, un marco espacial que no contiene ni datos ni información sobre CRS.

ext <- extent(raster)
crs(ext)

## [1] NA

class(ext)

## [1] "Extent"
## attr(,"package")
## [1] "raster"

4.2 Asignación de CRS a un objeto vectorial.

A continuación, transformaremos ese objeto extent, que sigue siendo un objeto ráster, en otro objeto espacial, pero ahora de la clase sf, que corresponde a un objeto vectorial. Para ello se utilizarán dos funciones:

• La función bbox convierte la “caja” extent en una caja similar bbox que puede ser leída por el paquetesf. Este objeto simplemente contiene las coordenadas del objeto.

caja <- st_bbox(ext)
caja

## xmin ymin xmax ymax 
##   -5   39   -2   42

• La función st_as_sfc transforma esas coordenadas en un objeto sf (un objeto vectorial, pero sólo sus dimensiones espaciales).

poligono <- st_as_sfc(caja)
rm(caja)
poligono

## Geometry set for 1 feature 
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -5 ymin: 39 xmax: -2 ymax: 42
## CRS:           NA

## POLYGON ((-5 39, -2 39, -2 42, -5 42, -5 39))

class(poligono)

## [1] "sfc_POLYGON" "sfc"

plot(poligono)

Sin embargo, este objeto carece de CRS. Para proporcionarle un CRS son posibles varias opciones

• Función st_crs(objeto) = "EPSG:4326"

st_crs(poligono) <- "EPSG:4326"
poligono

## Geometry set for 1 feature 
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -5 ymin: 39 xmax: -2 ymax: 42
## Geodetic CRS:  WGS 84

## POLYGON ((-5 39, -2 39, -2 42, -5 42, -5 39))

• Con la función st_set_crs(), que devuelve el objeto con el nuevo CRS:

borrame1 <- st_set_crs(poligono, 
                       "EPSG:4326")
borrame1

## Geometry set for 1 feature 
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -5 ymin: 39 xmax: -2 ymax: 42
## Geodetic CRS:  WGS 84

## POLYGON ((-5 39, -2 39, -2 42, -5 42, -5 39))

• Incluyendo el nuevo CRS como un argumento adicional en la función st_as_sf(... crs = "codigo EPSG").

borrame2 <- st_as_sf(poligono, 
                     "EPSG:4326")
borrame2

## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -5 ymin: 39 xmax: -2 ymax: 42
## Geodetic CRS:  WGS 84
##   "EPSG:4326"                              x
## 1   EPSG:4326 POLYGON ((-5 39, -2 39, -2 ...

rm(ext, borrame1, borrame2)

4.2.1 Reproyección del CRS de un objeto vectorial.

La función sf::st_transform(x, crs) transforma (“reproyecta”) cualquier objeto vectorial a un nuevo CRS.

Por ejemplo, al objeto polígono se le podría reasignar otro CRS (reproyección de CRS a CRS), por ejemplo la proyección MODIS Sinusoidal (EPSG: 6842; “+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs”.

poligono_sinuosoidal <- st_transform(poligono, 
                                     "EPSG:6842")

A continuación, representaremos ambos objetos al mismo tiempo.

# Una fila, dos columnas
par(mfcol = c(1, 2))
# Gráficos
plot(poligono, axes = TRUE)
plot(poligono_sinuosoidal, axes = TRUE)

# Volvemos al estado original
par(mfcol = c(1, 1))

Si quisiéramos superponer el objeto ráster y el objeto poligono habría que comprobar que poseen el mismo CRS. Si fuera así, habría que utilizar uno común. Recuérdese que raster no posee ningún CRS, mientras que acabamos de asignar un CRS a polígono.

utm <- crs("+init=epsg:23030")                   # Objeto con el CRS                
raster_ED50 <- projectRaster(raster,             # Reproyección del objeto raster (es necesario incluir argumento `crs =`)
                             crs = utm) 
poligono_ED50 <- st_transform(poligono,          # Reproyección del objeto vectorial
                              crs = utm)

Para verificar que la transformación se ha realizado correctamente, podemos superponer uno encima de otro.

plot(raster_ED50)
plot(poligono_ED50, add = TRUE)

Como no son necesarios para la siguiente sesión, se eliminarán los siguientes objetos del Global Environment

rm(list=ls())

5 EJEMPLO DESARROLLADO CON DATOS REALES

A continuación se revisan algunos procedimientos comunes para realizar estas transformaciones en R con datos reales. Primero, debe establecerse la carpeta dónde se encuentran los ficheros que serán utilizados en esta sesión.

setwd("D:/G174_OCW/Lab_3_preprocesamiento_imagenes_satelite/datos/proyeccion")
getwd()

## [1] "D:/G174_OCW/Lab_3_preprocesamiento_imagenes_satelite/datos/proyeccion"

Una vez creado dicha carpeta, se descargan los ficheros necesarios para esta actividad directamente desde el servidor web.

# download.file("https://personales.unican.es/rasillad/docencia/g174/preprocesamiento/datos/proyeccion.zip",     # Nombre del fichero a descargar
#              "./mi_carpeta/mi_fichero.zip",                                            # Ruta y nombre del fichero a descargar
#              method = "auto")

Y descomprimirse

#unzip("./mi_carpeta/mi_fichero.zip",                 # Ruta y nombre del fichero
#      exdir = "./mi_carpeta/")                      # Ruta de la carpeta que lo contiene

Un procedimiento alternativo es realizar todo el proceso manualmente. Para ello, los ficheros necesarios para esta práctica se descargan aquí

5.1 Objetos disponibles

En primer lugar se creará un fichero ráster A continuación se diseñará un objeto ráster formado por 81 celdas, con proyección geográfica.

mdt <- raster(ncol= 9,                         # Número de columnas del ráster
               nrow= 9,                        # Número de filas
               xmx= -3,xmn= -4.75,             # Coordenadas del eje X (considerémoslas la longitud)
               ymx= 41.5,ymn= 39.5)            # Coordenadas del eje y (idem latitud)
values(mdt) <- as.integer(runif(ncell(mdt),
                                 min = 1,
                                 max= 30))

A continuación se importará un fichero con los bordes administrativos de la Comunidad de Madrid.

comunidad <- st_read("D:/G174_OCW/Lab_3_preprocesamiento_imagenes_satelite/datos/proyeccion/comunidad.shp")

## Reading layer `comunidad' from data source 
##   `D:\G174_OCW\Lab_3_preprocesamiento_imagenes_satelite\datos\proyeccion\comunidad.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 365635.7 ymin: 4415272 xmax: 495483.9 ymax: 4557310
## Projected CRS: WGS 84 / UTM zone 30N

plot(comunidad, 
     col = "red")

También se importará el límite administrativo del municipio de Madrid.

library(rgeoboundaries)                            # Paquete con las unidades administrativas (todo el planeta)
library(dplyr)                                     # Paquete para filtrar/seleccionar unidades administrativas

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

municipio <- geoboundaries("Spain", "adm2")        # Solicitamos la extracción de los municipios ("adm") de España
municipio <- filter(municipio,                     # Filtramos el límite del municipio de Madrid
                  shapeName == "Madrid")
municipio <- municipio[ , -c(2:5)]                 # Eliminamos columnas innecesarias

Comprobamos el CRS de este nuevo objeto gráfico.

crs(municipio)

## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

Es habitual que en la mayoría de los paquetes que ofrecen datos espaciales el CRS básico corresponda a un sistema de coordenadas geográficas. La conversión a proyectadas es muy más fácil y genera menos errores.

Si se intenta representar la combinación del objeto ráster mdt, y superpuesto el objeto vectorial comunidad aparece lo siguiente.

plot(mdt) 
plot(st_geometry(comunidad), 
     col = "blue",
     add = TRUE)

Si se intenta realizar lo mismo con el objeto ráster y el objeto vectorial municipio el resultado es el siguiente.

plot(mdt) 
plot(st_geometry(municipio), 
     col = "red", 
     add = TRUE)

Por último, si se representan ambos límites administrativos, el resultado es el siguiente

plot(st_geometry(comunidad), 
     col = "lightgreen", 
     border = "gray") 
plot(st_geometry(municipio), 
     col = "red", add = TRUE)

El problema se debe a que esos objetos poseen diferentes CRS.

compareCRS(comunidad, municipio)

## [1] FALSE

compareCRS(mdt, municipio)

## [1] TRUE

Debemos conocer cuál es el CRS de cada objeto. En el caso del objeto ráster:

crs(mdt)

## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

En el caso de los objetos vectoriales

st_crs(comunidad)

## Coordinate Reference System:
##   User input: WGS 84 / UTM zone 30N 
##   wkt:
## PROJCRS["WGS 84 / UTM zone 30N",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 30N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-3,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     ID["EPSG",32630]]

st_crs(municipio)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

El objeto comunidad aparece en coordenadas proyectadas, mientras que el objeto municipio está proyectado en coordenadas geográficas, al igual que objeto ráster. La función plot()considera que las coordenadas del objeto comunidad se encuentran fuera del marco espacial en el que se inserta la capa municipio y, por tanto, no la representa, pero tampoco avisa del error ni proporciona un mensaje de advertencia. Esto también se hace evidente si comparamos:

• la extensión de cada objeto (le pediremos a R que muestra las coordenadas en formato vector para que no ocupe tanto):

as.vector (extent(comunidad))

## [1]  365635.7  495483.9 4415271.9 4557310.2

as.vector(extent(municipio))

## [1] -3.888957 -3.518235 40.312064 40.643333

as.vector(extent(mdt))

## [1] -4.75 -3.00 39.50 41.50

• Las proyecciones de ambos objetos:

projection(comunidad)

## [1] "+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs"

projection (municipio)

## [1] "+proj=longlat +datum=WGS84 +no_defs"

projection(mdt)

## [1] "+proj=longlat +datum=WGS84 +no_defs"

5.2 Transformación en un nuevo CRS (Reproyección)

Elegiremos como CRS común el sistema proyectado EPSG:23030 (ED50 / UTM zone 30N), que también es un sistema de coordenadas proyectadas. Hay varias posibilidades para realizar esta reproyección.

5.2.1 Reproyección de un objeto vectorial mediante la creación de un objeto CRS.

• Primero se puede crear un objeto con el nuevo sistema de proyección.

reproj <- "EPSG:23030"

A continuación, se transforma el sistema de referencia original del objeto municipio al nuevo CRS.

comunidad_ED50 <- st_transform(comunidad, 
                               crs = reproj)
crs(comunidad_ED50)

## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=30 +ellps=intl +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["ED50 / UTM zone 30N",
##     BASEGEOGCRS["ED50",
##         DATUM["European Datum 1950",
##             ELLIPSOID["International 1924",6378388,297,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4230]],
##     CONVERSION["UTM zone 30N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Europe - between 6°W and 0°W - Channel Islands (Jersey, Guernsey); France offshore; Gibraltar; Ireland offshore; Norway including Svalbard - offshore; Spain - onshore; United Kingdom - UKCS offshore."],
##         BBOX[35.26,-6,80.53,0]],
##     ID["EPSG",23030]]

Podríamos realizar lo mismo para reproyectar el objeto ráster.

mdt_ED50 <- projectRaster(mdt, 
                          crs = reproj)
crs(mdt_ED50)

## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=30 +ellps=intl +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["ED50 / UTM zone 30N",
##     BASEGEOGCRS["ED50",
##         DATUM["European Datum 1950",
##             ELLIPSOID["International 1924",6378388,297,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4230]],
##     CONVERSION["UTM zone 30N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Europe - between 6°W and 0°W - Channel Islands (Jersey, Guernsey); France offshore; Gibraltar; Ireland offshore; Norway including Svalbard - offshore; Spain - onshore; United Kingdom - UKCS offshore."],
##         BBOX[35.26,-6,80.53,0]],
##     ID["EPSG",23030]]

Ahora sí los límites administrativos coinciden.

plot(mdt_ED50) 
plot(st_geometry(comunidad), 
     add = TRUE)

5.2.2 Reproyección de un objeto vectorial con el CRS de un objeto ráster

Otra posibilidad es la reproyección de un objeto vectorial, en este caso, el objeto municipio con el CRS de un objeto raster, por ejemplo el objeto mdt_ED50.

El CRS del objeto vectorial es el siguiente.

crs(municipio)

## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

El CRS del objeto ráster aparece a continuación.

crs(mdt_ED50)

## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=30 +ellps=intl +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["ED50 / UTM zone 30N",
##     BASEGEOGCRS["ED50",
##         DATUM["European Datum 1950",
##             ELLIPSOID["International 1924",6378388,297,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4230]],
##     CONVERSION["UTM zone 30N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Europe - between 6°W and 0°W - Channel Islands (Jersey, Guernsey); France offshore; Gibraltar; Ireland offshore; Norway including Svalbard - offshore; Spain - onshore; United Kingdom - UKCS offshore."],
##         BBOX[35.26,-6,80.53,0]],
##     ID["EPSG",23030]]

A continuación, mediante la función sf::st_transform reproyectamos el objeto municipio al nuevo sistema de coordenadas, creando un nuevo objeto. En este caso, el argumento crs = crs(mdt_ED50) informa a R que tiene que utilizar como CRS el que corresponde al objeto mdt_ED50.

municipio_ED50 <- st_transform(municipio, 
                               crs = crs(mdt_ED50))

Ahora sí los límites administrativos coinciden.

plot(mdt_ED50) 
plot(st_geometry(comunidad_ED50), 
     col = "red", add = TRUE)
plot(st_geometry(municipio_ED50), 
     col = "yellow", add = TRUE)

Un último caso es la incorporación de otro objeto espacial que originalmente no dispone de proyección, porque ha sido importado como dataframe, por ejemplo, una serie de localidades (puntos). El procedimiento para asignar un CRS es algo más complicado.

Primero, se importa el fichero original (en formato .csv).

puntos <- read.csv("D:/G174_OCW/Lab_3_preprocesamiento_imagenes_satelite/datos/proyeccion/meteo_mad.csv")
str(puntos)

## 'data.frame':    9 obs. of  4 variables:
##  $ cdg     : int  82150 82190 82200 82210 82220 82230 82240 82260 82270
##  $ estacion: chr  " NAVACERRADA" " COLMENAR" " C. UNIVERSITARIA" " BARAJAS" ...
##  $ lon     : num  -4.02 -3.73 -3.72 -3.57 -3.68 ...
##  $ lat     : num  40.8 40.6 40.5 40.5 40.4 ...

class(puntos)

## [1] "data.frame"

Puede observarse que estos datos contienen información espacial al incluir dos columnas, una con la longitud de la localidad y otra con la latitud. Aprovechando estos datos, el dataframe se convertirá en un objeto espacial de la clase sf(un objeto vectorial). Dado que la localización de cada entidad está medida en grados, es necesario asignar como CRS el Sistema de Coordenadas Geográficas EPSG:4326 (WGS84).

puntos.sf <- st_as_sf(puntos, 
                      coords = c("lon", "lat"), 
                      crs = 4326)
class(puntos.sf)

## [1] "sf"         "data.frame"

A partir de aquí, podríamos transformar el CRS del objeto puntos.sf al sistema común.

meteo_ED50 <- st_transform(puntos.sf, 
                           crs = reproj)
meteo_ED50

## Simple feature collection with 9 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 414318.8 ymin: 4460848 xmax: 485993.5 ymax: 4517260
## Projected CRS: ED50 / UTM zone 30N
##     cdg          estacion                 geometry
## 1 82150       NAVACERRADA POINT (414318.8 4517260)
## 2 82190          COLMENAR POINT (438137.3 4500371)
## 3 82200  C. UNIVERSITARIA POINT (439309.5 4478159)
## 4 82210           BARAJAS POINT (452060.5 4482951)
## 5 82220            RETIRO POINT (442164.4 4474473)
## 6 82230    CUATRO VIENTOS POINT (433465.5 4469439)
## 7 82240            GETAFE POINT (438574.4 4460848)
## 8 82260       GUADALAJARA POINT (485993.5 4502013)
## 9 82270          TORREJON   POINT (462316 4483225)

Finalmente, se pueden cartografiar ambos conjuntos de datos.

plot(mdt_ED50) 
plot(st_geometry(comunidad_ED50), 
     col = "red", add = TRUE)
plot(st_geometry(municipio_ED50), 
     col = "yellow",
     border = "darkgreen", 
     add = TRUE) 
plot(st_geometry(meteo_ED50),
     col = "blue",
     add = TRUE)