Extract raster files from VIIRS VNP02 files provided in NetCDF format?

Question

I work in R with HDF files and sometimes i use HEG tool to deal with HDF files.

But now NASA has introduced something else that is NetCDF files with VIIRS VNP02 data. I want to extract I04 and I05 bands from a .nc file attached herewith. the extracted bands should be projected with UTM WGS84 Zone 45. Link to example NetCDF (.nc) file

Can anyone help me how to do this in R.

Answer 1

The spatial coordinates for the data in VNP02IMG.A2018228.1948.001.2018229023940.nc seem to appear in the VNP03IMG.A2018228.1948.001.2018229003536.nc file.

This can be downloaded via https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/5110/VNP03IMG/2018/228/VNP03IMG.A2018228.1948.001.2018229003536.nc but I'm not sure how you infer the last numbers from the VNP02 file. Note this matching is based on guesswork based on the two files having the same parts in the name and also basically the same geographic metadata and the same size grid:

$ gdalinfo VNP03IMG.A2018228.1948.001.2018229003536.nc | grep GRing
  GRingPointLatitude=28.95256 33.857861 54.62978 47.946575 
  GRingPointLongitude=103.91853 72.036797 72.617546 116.37637 
$ gdalinfo VNP02IMG.A2018228.1948.001.2018229023940.nc | grep GRing
  GRingPointLatitude=28.9526 33.857899 54.629799 47.946602 
  GRingPointLongitude=103.919 72.036797 72.6175 116.376 

Once got, you can read in the data using the HDF5 driver and point to the subdataset using the raster package:

I05 = raster(
 'HDF5:"VNP02IMG.A2018228.1948.001.2018229023940.nc"://observation_data/I05'
)

Then similar lets you read the coordinates into a raster grid:

lat = raster(
 'HDF5:"VNP03IMG.A2018228.1948.001.2018229003536.nc"://geolocation_data/latitude')
lon = raster(
 'HDF5:"VNP03IMG.A2018228.1948.001.2018229003536.nc"://geolocation_data/longitude')

These are a bit large for me to plot (6464 x 6400 pixels each). I'll take a sample:

> s = sample(prod(dim(I05)), 100000)
> d = data.frame(lon=lon[][s], lat=lat[][s], I05=I05[][s])
> ggplot(d,aes(x=lon, y=lat)) + geom_point(aes(col=I05))

Getting rid of the large values that flag various conditions makes a neater plot:

> d= d[d$I05 < 65500,]
> dim(d)
[1] 87104     3

There you can see the basic structure even though what you have is 100,000 random dots. The underlying "grid" structure is warped to that curvilinear outline.

If you want to work with this as a non-rectilinear grid you'll need the stars package and a lot of RAM to deal with the size. I'd also see if you can get an authoritative reference on this before doing any serious work with it.