Sistemas de Referencia de Coordenadas (SRC)

Dr. Francisco Zambrano

¿Qué hemos visto hasta el momento?

Manejo de tipos y estructuras de datos con
Datos Espaciales en
- Datos vectoriales ({sf})
- Datos raster ({terra})
Cargar y guardar datos espaciales
- Vectoriales
- Rasters

Contenidos clase

Sistemas de Referencias de Coordenadas (SRC)
- Datos Vectoriales
- Datos Raster

1. Manejo de CRS en R

Coordenadas geográficas: esféricas
Coordenadas proyectadas: cónica, cilíndrica y planas

Para manejar CRS en R se utiliza la librería PROJ

PROJ.4 fue usado por más de 30 años, luego mutó a \(PR \phi J\), dando paso a PROJ5, PROJ6,PROJ7,PROJ8. Ahora PROJ9 en 2022.

Así, hay tres formatos para manejar los CRS

proj.4
EPSG
Well-known Text (wkt)

{sf} y {terra} usan EPSG y wkt (well-known text)

2. Datos vectoriales: `{sf}`

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

2. Datos vectoriales: `{sf}`

Recomendación utilizar EPSG (European Petroleum Survey Group) y wkt para identificar el CRS, cuando sea posible. No se aconseja el uso de proj4.

Transformar de un CRS a otro con st_transform

nc = read_sf(system.file("shape/nc.shp", package="sf"))
head(nc$geometry[[4]][[1]][[1]],n=3)

          [,1]     [,2]
[1,] -76.00897 36.31960
[2,] -76.01735 36.33773
[3,] -76.03288 36.33598

nc.web_mercator <- st_transform(nc, 3857)
head(st_geometry(nc.web_mercator)[[4]][[1]][[1]],n=3)

         [,1]    [,2]
[1,] -8461242 4344697
[2,] -8462174 4347203
[3,] -8463903 4346960

3. Datos raster: `{terra}`

url <- 'https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_monthly/cogs/'
name <- 'chirps-v2.0.2021.09.cog'
(x <- rast(paste0('/vsicurl/',url,name)))

class       : SpatRaster 
dimensions  : 2000, 7200, 1  (nrow, ncol, nlyr)
resolution  : 0.05, 0.05  (x, y)
extent      : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : chirps-v2.0.2021.09.cog 
name        : chirps-v2.0.2021.09

3. Datos raster: `{terra}`

Para inspeccionar el SRC usamos la función crs

cat(crs(x))

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

3. Datos raster: `{terra}`

Asignar un CRS a un raster. Utilice crs

xx <- x
crs(xx) <- ''
crs(xx)

[1] ""

3. Datos raster: `{terra}`

Asignar un CRS a un raster. Utilice crs

xx2 <- xx
crs(xx) <- "+proj=longlat +datum=WGS84"
crs(xx2) <- "epsg:4326"

3. Datos raster: `{terra}`

Reproyectar

La reproyección de raster se pierde precisión.
Utiliza “nearest neighbor” para datos categóricos.
En otro caso utiliza interpolación “bilineal” por defecto.

r <- rast(xmin=-110, xmax=-90, ymin=40, ymax=60, ncols=40, nrows=40)
values(r) <- 1:ncell(r)
r

class       : SpatRaster 
dimensions  : 40, 40, 1  (nrow, ncol, nlyr)
resolution  : 0.5, 0.5  (x, y)
extent      : -110, -90, 40, 60  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
source(s)   : memory
name        : lyr.1 
min value   :     1 
max value   :  1600

3. Datos raster: `{terra}`

Reproyectar

plot(r)

3. Datos raster: `{terra}`

Reproyectar

newcrs <- "EPSG:32719"
pr1 <- terra::project(r, newcrs)
cat(crs(pr1))

PROJCRS["WGS 84 / UTM zone 19S",
    BASEGEOGCRS["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]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 19S",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-69,
            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",10000000,
            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["Navigation and medium accuracy spatial referencing."],
        AREA["Between 72°W and 66°W, southern hemisphere between 80°S and equator, onshore and offshore. Argentina. Bolivia. Brazil. Chile. Colombia. Peru."],
        BBOX[-80,-72,0,-66]],
    ID["EPSG",32719]]

Sistemas de Referencia de Coordenadas (SRC)

¿Qué hemos visto hasta el momento?

Contenidos clase

1. Manejo de CRS en R

2. Datos vectoriales: {sf}

2. Datos vectoriales: {sf}

3. Datos raster: {terra}

3. Datos raster: {terra}

3. Datos raster: {terra}

3. Datos raster: {terra}

3. Datos raster: {terra}

3. Datos raster: {terra}

3. Datos raster: {terra}

Recomendación

Generalmente tratar de transformar el SRC de un vectorial a un raster y no al revés.

2. Datos vectoriales: `{sf}`

2. Datos vectoriales: `{sf}`

3. Datos raster: `{terra}`

3. Datos raster: `{terra}`

3. Datos raster: `{terra}`

3. Datos raster: `{terra}`

3. Datos raster: `{terra}`

3. Datos raster: `{terra}`

3. Datos raster: `{terra}`

`Generalmente tratar de transformar el SRC de un vectorial a un raster y no al revés.`