Some problems with the netcdf-cf mapper

[mortenwh]

Demo file: https://thredds.met.no/thredds/fileServer/aromearcticarchive/2019/06/15/arome_arctic_vtk_20190615T15Z.nc (download with, e.g., wget)

• The projection based on a subdataset is slightly different from the one based on the projection variable (projection_lambert) provided in the netcdf file. It is probably best to use the projection variable if it is present.

• During work with the above file, I discovered that the selection of bands based on the netcdf_dim kwarg is wrong. The time selection works, but it is not possible to filter on the pressure dimension (ref above file).

• Nansat fails to open a downloaded S1A file (wget http://nbstds.met.no/thredds/fileServer/NBS/S1A/2019/06/15/IW/S1A_IW_GRDM_1SDV_20190615T151935_20190615T152004_027692_03202F_6E58.nc). Also, the scatterometer mapper is selected by Nansat, although it should be the netcdf-cf mapper. We may need to make a specific S1 netcdf-cf mapper.

• If I reproject the arome dataset to an S1 dataset (http://nbstds.met.no/thredds/dodsC/NBS/S1A/2019/06/15/IW/S1A_IW_GRDM_1SDV_20190615T151935_20190615T152004_027692_03202F_6E58.nc) available from opendap), the bands cannot be read. Changing the projection as suggested above (point 1) makes it readable, but all the data is NaN. It works fine if I use the opendap version of the arome dataset (https://thredds.met.no/thredds/dodsC/aromearcticarchive/2019/06/15/arome_arctic_vtk_20190615T15Z.nc).

• Also check S2 netcdf datasets

[mortenwh]

Regarding reading the netcdf's from satellittdata.no:

OPeNDAP (http://nbstds.met.no/thredds/dodsC/NBS/S1A/2019/06/15/IW/S1A_IW_GRDM_1SDV_20190615T151935_20190615T152004_027692_03202F_6E58.nc):

unavailable at the moment of writing the comment - will edit later...

Downloaded NetCDF:

In [1]: import gdal                                                                                                              

In [2]: import time                                                                                                              

In [3]: gds = gdal.Open('NETCDF:"S1A_IW_GRDM_1SDV_20190615T151935_20190615T152004_027692_03202F_6E58.nc":Amplitude_VV')          

In [4]: start, amp0, end = time.time(), gds.ReadAsArray(), time.time()                                                           

In [5]: end-start # seconds                                                                                                      
Out[5]: 598.0722901821136

In [6]: (end-start)/60 # minutes                                                                                                 
Out[6]: 9.967871503035228

In [7]: from netCDF4 import Dataset                                                                                              

In [8]: ds = Dataset('S1A_IW_GRDM_1SDV_20190615T151935_20190615T152004_027692_03202F_6E58.nc')                                   

In [9]: start, amp, end = time.time(), ds.variables['Amplitude_VV'][:].data, time.time()                                         

In [10]: end-start # seconds                                                                                                     
Out[10]: 0.4655623435974121

10 minutes with gdal vs 0.5 seconds with netCDF4!

Then, transform the netcdf file:

gdal_translate -of netcdf -co "WRITE_BOTTOMUP=NO" -sds S1A_IW_GRDM_1SDV_20190615T151935_20190615T152004_027692_03202F_6E58.nc out.nc
nccopy -7 -c x/3,y/3 out_04.nc out_04_chunked.nc

This shortens the gdal reading time (but now there is only one band in out*.nc):

In [7]: gds = gdal.Open('out_04_chunked.nc')                                                                                     

In [8]: start, amp, end = time.time(), gds.ReadAsArray(), time.time()                                                            

In [9]: end-start # seconds                                                                                                      
Out[9]: 3.566657304763794

In [10]:                                                                                                                         

In [10]: gds = gdal.Open('out_04.nc')                                                                                            

In [11]: start, amp, end = time.time(), gds.ReadAsArray(), time.time()                                                           

In [12]: end-start # seconds                                                                                                     
Out[12]: 2.888197183609009

In [13]: 

ln [14]: ds = Dataset('out_04.nc')                                                                                               

In [15]: start, amp, end = time.time(), ds.variables['Amplitude_VV'][:].data, time.time()                                        

In [16]: end-start # seconds                                                                                                     
Out[16]: 1.4296467304229736

netCDF4 takes more time to read the new file compared to the original but is still much faster than gdal...

[mortenwh]

Reprojection of the Arome dataset now works with the new updates using the downloaded netCDF files:

In [1]: from nansat.nansat import Nansat                                                                                         

In [2]: import numpy as np                                                                                                       

In [3]: a0fn='arome_arctic_vtk_20190615T15Z.nc'                                                                                  

In [4]: fn = 'S1A_IW_GRDM_1SDV_20190615T151935_20190615T152004_027692_03202F_6E58.nc'                                            

In [5]: a0n = Nansat(a0fn, bands=['x_wind', 'y_wind'], netcdf_dim={'time': np.datetime64('2019-06-15T15:00')})                   

In [6]: n = Nansat(fn)                                                                                                           
=>Sentinel-1 GRD data<=
/home/vagrant/Miniconda3-4.6.14-Linux-x86_64/envs/py3openwind/lib/python3.6/site-packages/nansat/mappers/mapper_netcdf_cf.py:365: UserWarning: GDAL cannot determine the dataset projection - using Nansat spatial reference WKT, assuming a regular longitude/latitude grid
  warnings.warn('GDAL cannot determine the dataset projection - using Nansat ' \
/home/vagrant/Miniconda3-4.6.14-Linux-x86_64/envs/py3openwind/lib/python3.6/site-packages/nansat/mappers/mapper_netcdf_cf.py:304: UserWarning: KeyError: 'NETCDF_DIM_time' - Mapping key not found. Continuing without time metadata for band lat
  band_metadata['NETCDF_VARNAME']))
/home/vagrant/Miniconda3-4.6.14-Linux-x86_64/envs/py3openwind/lib/python3.6/site-packages/nansat/mappers/mapper_netcdf_cf.py:304: UserWarning: KeyError: 'NETCDF_DIM_time' - Mapping key not found. Continuing without time metadata for band lon
  band_metadata['NETCDF_VARNAME']))
/home/vagrant/Miniconda3-4.6.14-Linux-x86_64/envs/py3openwind/lib/python3.6/site-packages/nansat/mappers/mapper_netcdf_cf.py:304: UserWarning: KeyError: 'NETCDF_DIM_time' - Mapping key not found. Continuing without time metadata for band swathList
  band_metadata['NETCDF_VARNAME']))

In [7]: a0n.reproject(n)                                                                                                         

In [8]: a0n['x_wind_10m']                                                                                                        
Out[8]: 
array([[ 9.09290028,  9.09290028,  9.09290028, ...,  9.08801746,
         9.08801746,  9.08801746],
       [ 9.09290028,  9.09290028,  9.09290028, ...,  9.08801746,
         9.08801746,  9.08801746],
       [ 9.09290028,  9.09290028,  9.09290028, ...,  9.08801746,
         9.08801746,  9.08801746],
       ..., 
       [ 5.96741199,  5.96741199,  5.96741199, ...,  8.51916981,
         8.51916981,  8.51916981],
       [ 5.96741199,  5.96741199,  5.96741199, ...,  8.51916981,
         8.51916981,  8.51916981],
       [ 5.96741199,  5.96741199,  5.96741199, ...,  8.51916981,
         8.51916981,  8.51916981]], dtype=float32)

In [9]: n['Amplitude_VV']
Out[9]: 
array([[54, 49, 50, ...,  0,  0,  0],
       [67, 60, 58, ...,  0,  0,  0],
       [75, 68, 67, ...,  0,  0,  0],
       ..., 
       [ 0,  0,  0, ...,  0,  0,  0],
       [ 0,  0,  0, ...,  0,  0,  0],
       [ 0,  0,  0, ...,  0,  0,  0]], dtype=uint16)

[mortenwh]

S2 file: wget http://nbstds.met.no/thredds/fileServer/NBS/S2A/2019/07/01/S2A_MSIL1C_20190701T160911_N0207_R140_T29XNK_20190701T180203.nc