Manera eficiente de trazar datos en una cuadrícula irregular

Question

Feb 23, 2018, 05:36 PM

Manera eficiente de trazar datos en una cuadrícula irregular

Trabajo con datos satelitales organizados en una cuadrícula bidimensional irregular cuyas dimensiones son la línea de exploración (a lo largo de la dimensión de la pista) y el píxel del suelo (a través de la dimensión de la pista). La información de latitud y longitud para cada píxel central se almacena en variables de coordenadas auxiliares, así como los pares de coordenadas de las cuatro esquinas (las coordenadas de latitud y longitud se dan en el elipsoide de referencia WGS84). Los datos se almacenan en archivos netCDF4.

Lo que estoy tratando de hacer es trazar eficientemente estos archivos (y posiblemente una combinación de archivos, ¡el siguiente paso!) En un mapa proyectado.

Mi enfoque hasta ahora, inspirado porJeremy Voisey respuesta a estopregunta, ha sido crear un marco de datos que vincule mi variable de interés a los límites de píxeles y usarggplot2 congeom_polygon para la trama real.

Permítanme ilustrar mi flujo de trabajo y me disculpo de antemano por el enfoque ingenuo: acabo de comenzar a codificar con R desde hace una o dos semanas.

Nota

Para reproducir completamente el problema:
1. descargue los dos marcos de datos:so2df.Rda (22M) ypixel_corners.Rda (26M)
2. cárguelos en su entorno, p.

so2df <- readRDS(file="so2df.Rda")
pixel_corners <- readRDS(file="pixel_corners.Rda")

salte al paso "Fusionar los marcos de datos".Configuración inicial

Voy a leer los datos y los límites de latitud / longitud de mi archivo.

library(ncdf4)
library(ggplot2)
library(ggmap) 
# set path and filename
ncpath <- "/Users/stefano/src/s5p/products/e1dataset/L2__SO2/"
ncname <- "S5P_OFFL_L2__SO2____20171128T234133_20171129T003956_00661_01_022943_00000000T000000"  
ncfname <- paste(ncpath, ncname, ".nc", sep="")
nc <- nc_open(ncfname)

# save fill value and multiplication factors
mfactor = ncatt_get(nc, "PRODUCT/sulfurdioxide_total_vertical_column", 
                    "multiplication_factor_to_convert_to_DU")
fillvalue = ncatt_get(nc, "PRODUCT/sulfurdioxide_total_vertical_column", 
                      "_FillValue")

# read the SO2 total column variable
so2tc <- ncvar_get(nc, "PRODUCT/sulfurdioxide_total_vertical_column")

# read lat/lon of centre pixels
lat <- ncvar_get(nc, "PRODUCT/latitude")
lon <- ncvar_get(nc, "PRODUCT/longitude")

# read latitude and longitude bounds
lat_bounds <- ncvar_get(nc, "GEOLOCATIONS/latitude_bounds")
lon_bounds <- ncvar_get(nc, "GEOLOCATIONS/longitude_bounds")

# close the file
nc_close(nc)
dim(so2tc)
## [1]  450 3244

Entonces, para este archivo / pase, hay 450 píxeles de tierra para cada una de las 3244 líneas de exploración.

Crear los marcos de datos

Aquí creo dos marcos de datos, uno para los valores, con algo de procesamiento posterior, y uno para los límites de latitud / longitud, luego combino los dos marcos de datos.

so2df <- data.frame(lat=as.vector(lat), lon=as.vector(lon), so2tc=as.vector(so2tc))
# add id for each pixel
so2df$id <- row.names(so2df)
# convert to DU
so2df$so2tc <- so2df$so2tc*as.numeric(mfactor$value)
# replace fill values with NA
so2df$so2tc[so2df$so2tc == fillvalue] <- NA
saveRDS(so2df, file="so2df.Rda")
summary(so2df)

##       lat              lon              so2tc              id           
##  Min.   :-89.97   Min.   :-180.00   Min.   :-821.33   Length:1459800    
##  1st Qu.:-62.29   1st Qu.:-163.30   1st Qu.:  -0.48   Class :character  
##  Median :-19.86   Median :-150.46   Median :  -0.08   Mode  :character  
##  Mean   :-13.87   Mean   : -90.72   Mean   :  -1.43                     
##  3rd Qu.: 31.26   3rd Qu.: -27.06   3rd Qu.:   0.26                     
##  Max.   : 83.37   Max.   : 180.00   Max.   :3015.55                     
##                                     NA's   :200864

Guardé este marco de datos comoso2df.Rda aquí (22M).

num_points = dim(lat_bounds)[1]
pixel_corners <- data.frame(lat_bounds=as.vector(lat_bounds), lon_bounds=as.vector(lon_bounds))
# create id column by replicating pixel's id for each of the 4 corner points
pixel_corners$id <- rep(so2df$id, each=num_points)
saveRDS(pixel_corners, file="pixel_corners.Rda")
summary(pixel_corners)


##    lat_bounds       lon_bounds           id           
##  Min.   :-89.96   Min.   :-180.00   Length:5839200    
##  1st Qu.:-62.29   1st Qu.:-163.30   Class :character  
##  Median :-19.86   Median :-150.46   Mode  :character  
##  Mean   :-13.87   Mean   : -90.72                     
##  3rd Qu.: 31.26   3rd Qu.: -27.06                     
##  Max.   : 83.40   Max.   : 180.00

Como se esperaba, el marco de datos de límites lat / lon es cuatro veces más grande que el marco de datos de valor (cuatro puntos por cada píxel / valor).
Guardé este marco de datos comopixel_corners.Rda aquí (26M).

Fusionar los marcos de datos

Luego combino los dos marcos de datos por id:

start_time <- Sys.time()
so2df <- merge(pixel_corners, so2df, by=c("id"))
time_taken <- Sys.time() - start_time
print(paste(dim(so2df)[1], "rows merged in", time_taken, "seconds"))

## [1] "5839200 rows merged in 42.4763631820679 seconds"

Como puede ver, es un proceso intensivo de CPU de alguna manera. Me pregunto qué pasaría si trabajara con 15 archivos a la vez (cobertura global).

Trazando los datos

Ahora que tengo mis esquinas de píxeles vinculadas al valor de píxeles, puedo trazarlas fácilmente. Por lo general, estoy interesado en una región particular de la órbita, por lo que hice una función que subconjusta el marco de datos de entrada antes de trazarlo:

PlotRegion <- function(so2df, latlon, title) {
  # Plot the given dataset over a geographic region.
  #
  # Args:
  #   df: The dataset, should include the no2tc, lat, lon columns
  #   latlon: A vector of four values identifying the botton-left and top-right corners 
  #           c(latmin, latmax, lonmin, lonmax)
  #   title: The plot title

  # subset the data frame first
  df_sub <- subset(so2df, lat>latlon[1] & lat<latlon[2] & lon>latlon[3] & lon<latlon[4])

  subtitle = paste("#Pixel =", dim(df_sub)[1], "- Data min =", 
                   formatC(min(df_sub$so2tc, na.rm=T), format="e", digits=2), "max =", 
                   formatC(max(df_sub$so2tc, na.rm=T), format="e", digits=2))

  ggplot(df_sub) + 
    geom_polygon(aes(y=lat_bounds, x=lon_bounds, fill=so2tc, group=id), alpha=0.8) +
    borders('world', xlim=range(df_sub$lon), ylim=range(df_sub$lat), 
            colour='gray20', size=.2) + 
    theme_light() + 
    theme(panel.ontop=TRUE, panel.background=element_blank()) +
    scale_fill_distiller(palette='Spectral') +
    coord_quickmap(xlim=c(latlon[3], latlon[4]), ylim=c(latlon[1], latlon[2])) +
    labs(title=title, subtitle=subtitle, 
         x="Longitude", y="Latitude", 
         fill=expression(DU)) 
}

Luego invoco mi función sobre una región de interés, por ejemplo, veamos qué sucede en Hawai:

latlon = c(17.5, 22.5, -160, -154)
PlotRegion(so2df, latlon, expression(SO[2]~total~vertical~column))

Ahí están, mis píxeles, y lo que parece ser un penacho de SO2 del Mauna Loa. Por favor ignore los valores negativos por ahora. Como puede ver, el área de los píxeles varía hacia el borde de la franja (diferente esquema de agrupamiento).

Traté de mostrar la misma trama sobre los mapas de Google, usando ggmap:

PlotRegionMap <- function(so2df, latlon, title) {
  # Plot the given dataset over a geographic region.
  #
  # Args:
  #   df: The dataset, should include the no2tc, lat, lon columns
  #   latlon: A vector of four values identifying the botton-left and top-right corners 
  #           c(latmin, latmax, lonmin, lonmax)
  #   title: The plot title

  # subset the data frame first
  df_sub <- subset(so2df, lat>latlon[1] & lat<latlon[2] & lon>latlon[3] & lon<latlon[4])

  subtitle = paste("#Pixel =", dim(df_sub)[1], "Data min =", formatC(min(df_sub$so2tc, na.rm=T), format="e", digits=2), 
                   "max =", formatC(max(df_sub$so2tc, na.rm=T), format="e", digits=2))
  base_map <- get_map(location = c(lon = (latlon[4]+latlon[3])/2, lat = (latlon[1]+latlon[2])/2), zoom = 7, maptype="terrain", color="bw")

  ggmap(base_map, extent = "normal")  +
    geom_polygon(data=df_sub, aes(y=lat_bounds, x=lon_bounds,fill=so2tc, group=id),  alpha=0.5) +
    theme_light() + 
    theme(panel.ontop=TRUE, panel.background=element_blank()) +
    scale_fill_distiller(palette='Spectral') +
    coord_quickmap(xlim=c(latlon[3], latlon[4]), ylim=c(latlon[1], latlon[2])) +
    labs(title=title, subtitle=subtitle, 
         x="Longitude", y="Latitude", 
         fill=expression(DU)) 

}

Y esto es lo que obtengo:

latlon = c(17.5, 22.5, -160, -154)
PlotRegionMap(so2df, latlon, expression(SO[2]~total~vertical~column))

Preguntas¿Hay alguna forma más eficiente de abordar este problema? Estaba leyendo sobre elsf paquete, y me preguntaba si podría definir un marco de datos de puntos (valores + coordenadas de píxel central), y tenersf inferir automáticamente los límites de píxeles. Eso me salvaría de tener que confiar en los límites lat / lon definidos en mi conjunto de datos original y tener que fusionarlos con mis valores. Podría aceptar una pérdida de precisión en las áreas de transición hacia el borde de la franja, de lo contrario la cuadrícula es bastante regular, cada píxel es de 3.5x7 km ^ 2 de gran tamaño.¿Volver a cuadricular mis datos en una cuadrícula normal (cómo?), Posiblemente agregando píxeles adyacentes, mejoraría el rendimiento? Estaba pensando en usar elraster paquete, que, como entendí, requiere datos en una cuadrícula regular. Esto debería ser útil a escala global (por ejemplo, diagramas en Europa), donde no necesito trazar los píxeles individuales; de hecho, ni siquiera podría verlos.¿Necesito volver a proyectar mis datos al trazar sobre Google Map?

[preguntas cosméticas adicionales]

¿Existe una forma más elegante de subconjugar mi marco de datos en una región identificada por sus cuatro puntos de esquina?¿Cómo podría cambiar la escala de colores para que los valores más altos se destaquen con respecto a los valores más bajos? He experimentado con la escala logarítmica con malos resultados.