R: ¿Debo transformar mi Raster en objetos espaciales aquí?

Question

Esta pregunta es un seguimiento basado en algunos de los consejos proporcionados en las respuestas y comentarios en Conversión de ráster a puntos espaciales sin pérdida de información en R .

En primer lugar, para proporcionar un poco de contexto, tengo un área geográfica de interés (NSW en Australia), y también tengo un shapefile que contiene las ubicaciones de los drenajes en la misma área. Tengo un ráster que contiene píxeles en la zona que me interesa (este ráster se almacena en una variable llamada ndvi_median ).

Mi objetivo es crear un ráster en el que los píxeles sean el área de interés, y en el que el atributo del ráster sea la distancia más corta (en m, aunque los km también están bien) entre cada píxel y el desagüe más cercano. Afortunadamente, ya hubo una pregunta similar en este intercambio de pilas: Encontrar la distancia entre los píxeles de la trama y las características de la línea en R .

La única diferencia es que mis shapefiles y rasters están en lat-long, pero yo quiero las distancias en metros. He leído en otro post de stack exchange que hay que proyectar a otro CRS - voy a usar +proj=lcc +lat_1=-30.75 +lat_2=-35.75 +lat_0=-33.25 +lon_0=147 +x_0=9300000 +y_0=4500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs para esta tarea. Por lo tanto, mi enfoque general es:

1. Reproyectar cuando sea necesario utilizando un nuevo SIR:

   ndvi_median # contains pixels in the geographic area of interest class : RasterLayer dimensions : 1122, 181, 203082 (nrow, ncol, ncell) resolution : 0.09978817, 0.008333333 (x, y) extent : 140.9993, 159.0609, -37.50702, -28.15702 (xmin, xmax, ymin, ymax) crs : +proj=longlat +ellps=GRS80 +no_defs source : C:/Users/thoma/AppData/Local/Temp/RtmpWGYBpT/raster/r_tmp_2021-05-18_120505_15536_85040.grd names : layer values : -1966, 9009.5 (min, max)

   GeoDat_NSW ## Contains information on drainages class : SpatialLinesDataFrame features : 1868 extent : 140.9993, 153.5635, -37.50702, -28.15702 (xmin, xmax, ymin, ymax) crs : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +nodefs variables : 12 names : FNODE, TNODE, LPOLY, RPOLY, LENGTH, AUS25DGD, ID, FEAT_CODE, NAME, PERENNIAL, Q_INFO, UFI min values : 5952, 5952, 0, 0, 0.003601388621079, 6150, 6150, canal, ABERCROMBIE RIVER, 0, BI000001, DJ00010416 max values : 9039, 8972, 0, 0, 3.84625259743069, 8737, 8737, watercours_l, YOWRIE RIVER, 2, BJ000004, DJ00013003

   Spat_points_NSW<-projectRaster(ndvi_median,crs="+proj=lcc +lat_1=-30.75 +lat_2=-35.75 +lat_0=-33.25 +lon_0=147 +x_0=9300000 +y_0=4500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs") Spat_points_NSW_2<-as(Spat_points_NSW,"SpatialPoints") GeoDat_Line<-as(GeoDat_NSW,"SpatialLines") GeoDat_Line_transformed<-spTransform(GeoDat_Line,CRS("+proj=lcc +lat_1=-30.75 +lat_2=-35.75 +lat_0=-33.25 +lon_0=147 +x_0=9300000 +y_0=4500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))

2. Ahora que todo ha sido reproyectado utilizando el nuevo SIR, aplico el método sugerido en Encontrar la distancia entre los píxeles de la trama y las características de la línea en R para encontrar la distancia entre las celdas y la línea de drenaje más cercana. De hecho, esta es una razón clave por la que sentí la necesidad de convertir todo en objetos sp (que se requieren en los argumentos del gdistance función).

   Dist_drainage<-gDistance(Spat_points_NSW_2,GeoDat_Line_transformed,byid=TRUE) # finds distance between each cell and every drainage line Min_Dist<-apply(Dist_drainage,2,min) # finds the minimum distance for each cell Spat_points_NSW_3<-Spat_points_NSW Spat_points_NSW_3[]=Min_Dist # Effectively creating a raster with the distance information as the raster attributes

3. Ahora, sólo tengo que convertir las coordenadas en Spat_points_NSW_3 de nuevo en lat-long.

   Reproject to lat-long

   Spat_points_NSW_4<-projectRaster(Spat_points_NSW_3,crs=crs(ndvi_median))

Desgraciadamente, este método no funciona del todo porque en la línea de abajo me sale un error

> Min_Dist<-apply(Dist_drainage,2,min)
Error: cannot allocate vector of size 1.7 Gb

Esto me lleva a preguntarme si existe una forma alternativa de abordar este problema. He visto una solución prometedora utilizando geosphere::dist2Line en Calcular la distancia entre los puntos y el polígono más cercano en R Pero, lamentablemente, el método tarda demasiado en ejecutarse (aunque sólo lo detuve después de ejecutarlo durante unos 40 minutos; puedo intentar dejarlo en funcionamiento toda la noche si me parece el método más óptimo).

Answer 1

Sería de gran ayuda que proporcionaras un ejemplo mínimo, autocontenido y reproducible.

Este es un enfoque con terra que funciona con los datos del ejemplo. No lo he probado con conjuntos de datos más grandes. Pero rasteriza las celdas tocadas por la línea, y luego calcula la distancia al centro de estas celdas; así que introduce algún error.

library(terra)
r <- rast(ncol=36,nrow=18, crs="+proj=longlat +datum=WGS84")
v <- vect(matrix(c(-120, 47, -137, 9, -124, -12, -129, -41), ncol=2, byrow=TRUE), type="lines", crs=crs(r))

d <- distance(r, v) 
plot(d / 1000)
lines(v)

En contra de la creencia popular, la transformación de los datos lon/lat a un SIR plano no es un buen enfoque para medir la distancia (o el área), a menos que se utilice un SIR plano que no distorsione la distancia (o que sea de igual área). Pero es cierto que geos::gDistance lo requiere. Implementaré el método de la geosfera en terra para que se ejecute (espero que mucho) más rápido.