Superponer raster (archivos h5) con archivos shape para obtener datos de luz nocturna en R

Question

Estoy tratando de combinar ".h5" archivos de datos de luz nocturna de Hawai con zipcode-shapefile de Hawai, pero estoy recibiendo el siguiente error:

#Error: [crop] extents do not overlap

Los códigos R están a continuación, y los datos se pueden descargar desde google drive con el enlace de abajo. https://drive.google.com/drive/u/2/folders/1IpkbgITblqxU9sArUFtD_-Oo3WtF6wQg

Nota: Si prefiere utilizar el enlace oficial para los datos (enlace más abajo), deberá registrarse.

library(terra) # the latest version of raster
library(sf)
# get hc files for jan 1 2022 of hawaii from the link below (need to register for getting data):
# https://ladsweb.modaps.eosdis.nasa.gov/search/order/4/VNP46A2--5000/2022-01-01..2022-01-01/N/-159.9,22.4,-154.4,18.8
# files are VNP46A2.A2022001.h02v07.001.2022009102111.h5 and 
# VNP46A2.A2022001.h02v06.001.2022009101742.h5
# read two hc tile files for day 1 of 2022
temp_files<-list.files(pattern="VNP46A2.A2022001")

# read each hc files with rast
temp_files_rast <- sprc(lapply(temp_files, rast))

# merge raster files
data_raster <- mosaic(temp_files_rast)

# we need `DNB_BRDF-Corrected_NTL` layer only
data_raster_req <- data_raster$`DNB_BRDF-Corrected_NTL`

# assign crf since no crf is provided
crs(data_raster_req) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" 

# we need to get the nightlight data for zipcode level, so get the shape file to overlay on raster file
# https://files.hawaii.gov/dbedt/op/gis/data/zcta20.shp.zip

haw_sf<-read_sf("zcta20.shp") 

# convert shape file to spatial class object
haw_sp<-as(haw_sf,"Spatial") # gives spatial class object

# to make sure shape file overlays with raster
haw_sp<-sp::spTransform(haw_sp,crs(data_raster_req))

# convert sp to sf
haw_sf<-as(haw_sp,"sf")

# next crop raster file for area of my interest
data_raster_crop <- crop(data_raster_req, haw_sp)

Answer 1

Cuando lees los archivos HDF, no obtienes ninguna coordenada de latitud y longitud, y hay mensajes de advertencia que te lo indican y que pareces haber ignorado:

temp_files_rast <- sprc(lapply(temp_files, rast))
Warning messages:
1: [rast] unknown extent

2: [rast] unknown extent

Si cargas una trama podrás ver un poco más:

> r1 = rast("./VNP46A2.A2022001.h02v06.001.2022009101742.h5")
Warning message:
[rast] unknown extent

> r1
class       : SpatRaster 
dimensions  : 2400, 2400, 7  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0, 2400, 0, 2400  (xmin, xmax, ymin, ymax)

Así pues, a los rásteres se les asigna una extensión en píxeles (se trata de rásteres de 2400x2400 píxeles).

Cuando se le asigna un CRS lat-long (no "crf") se está diciendo que este raster cubre el área de 0 latitud a 2400 latitud, lo que no tiene sentido.

A partir de ahí todo se tuerce. Realmente necesitas saber la extensión correcta en lat-long de los archivos HDF, y no veo esto en los metadatos HDF.

Si supone que el HDF y el shapefile representan la misma extensión, puede configurarlos para que tengan la misma extensión y CRS:

ext(data_raster_req) = ext(haw_sf)
crs(data_raster_req) = crs(haw_sf)

que produce esto cuando se mapea:

> plot(data_raster_req)
> plot(st_geometry(haw_sf),add=TRUE)

lo que podría ser correcto - no lo sé. Pero así es como se hace.

Pero esencialmente, necesita la extensión de los archivos HDF.

Answer 2

Como ha mencionado @Spacedman, sus datos ráster carecen de una definición adecuada del sistema de referencia de coordenadas y de la extensión.

Según la información sobre el producto para VNP46A2, su producto tiene una cobertura global y la resolución espacial se especifica con 15 segundos de arco. Los datos vienen en mosaicos (aquí: h02v06 y h02v07), así que todo lo que necesitas es encontrar la extensión de los mosaicos utilizados. He encontrado un shapefile con las geometrías de los mosaicos, por ejemplo, en blackmarble.gsfc.nasa.gov .

Utilización de {terra} estas propiedades se pueden (extraer y) definir así:

library(terra)
#> terra 1.6.17

# read hdf5 file
r <- rast("VNP46A2.A2022001.h02v06.001.2022009101742.h5")
#> Warning: [rast] unknown extent
r
#> class       : SpatRaster 
#> dimensions  : 2400, 2400, 7  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 2400, 0, 2400  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> sources     : VNP46A2.A2022001.h02v06.001.2022009101742.h5://DNB_BRDF-Corrected_NTL  
#>               VNP46A2.A2022001.h02v06.001.2022009101742.h5://DNB_Lunar_Irradiance  
#>               VNP46A2.A2022001.h02v06.001.2022009101742.h5://Gap_Filled_DNB_BRDF-Corrected_NTL  
#>               ... and 4 more source(s)
#> varnames    : DNB_BRDF-Corrected_NTL 
#>               DNB_Lunar_Irradiance 
#>               Gap_Filled_DNB_BRDF-Corrected_NTL 
#>               ...
#> names       : DNB_B~d_NTL, DNB_L~iance, Gap_F~d_NTL, Lates~ieval, Manda~_Flag, QF_Cl~_Mask, ... 
#> min values  :           1,          ? ,          ? ,          ? ,          ? ,          ? , ... 
#> max values  :       65535,          ? ,          ? ,          ? ,          ? ,          ? , ...

# read tile grid
tiles <- vect("BlackMarbleTiles.shp")

# get relevant index
ind <- tiles[["TileID"]] == "h02v06"

# subset grid to one tile, get extent
e <- tiles[ind] |> ext()

# adjust crs and extent of the tile in use
crs(r) <- "epsg:4326"
ext(r) <- e

# subset data to relevant layer and inspect
r[["DNB_BRDF-Corrected_NTL"]]
#> class       : SpatRaster 
#> dimensions  : 2400, 2400, 1  (nrow, ncol, nlyr)
#> resolution  : 0.004166667, 0.004166667  (x, y)
#> extent      : -160, -150, 20.00012, 30.00012  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : VNP46A2.A2022001.h02v06.001.2022009101742.h5://DNB_BRDF-Corrected_NTL 
#> varname     : DNB_BRDF-Corrected_NTL 
#> name        : DNB_BRDF-Corrected_NTL 
#> min value   :                      1 
#> max value   :                  65535

# resolution equal to 15 arc-seconds?
res(r)[1] == 1 / 3600 * 15
#> [1] TRUE