3

I am trying to combine ".h5" nightlight data files of Hawaii with zipcode-shapefile from Hawaii, but I am getting the following error:

#Error: [crop] extents do not overlap

The R codes are below, and the data can be downloaded from the google drive with the link below. https://drive.google.com/drive/u/2/folders/1IpkbgITblqxU9sArUFtD_-Oo3WtF6wQg

Note: If you prefer to use the the official link for data (link below), you need to register.

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)

2 Answers 2

2

When you read the HDF files, you aren't getting any coordinates for latitude and longitude, and there are warning messages telling you this which you seem to have ignored:

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

If you load one raster you can see a bit more:

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

So the rasters are given an extent in pixels (these are 2400x2400 pixel rasters).

When you then assign it a lat-long CRS (not "crf") you are saying this raster covers the area from 0 latitude to 2400 latitude, which is meaningless.

From then on it all goes wrong. You really need to know the correct extent in lat-long of the HDF files, and I don't see this in the HDF metadata.

If you do assume the HDF and the shapefile represent the same extent you can set them to have the same extent and CRS:

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

which produces this when mapped:

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

enter image description here

which might be right - I don't know. But that's how you do it.

But essentially, you need the extent of the HDF files.

2

As @Spacedman has mentioned, your raster data is lacking a proper definition of the coordinate reference system and extent.

According to the product information for VNP46A2, your product has a global coverage and the spatial resolution is specified with 15 arc-seconds. Data comes in tiles (here: h02v06 and h02v07), so all you need is to find the extent of the tiles used. I found a shapefile with tile geometries e.g. at blackmarble.gsfc.nasa.gov.

Making use of {terra}, these properties can be (extracted and) defined like this:

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.