python-netcdf4: Writing file leaves CRS as 'undefined' and extent not correct

Question

I am trying to write a netcdf file in my file system and then open it in QGIS for visualization purposes. I have 3 numpy arrays with size 234 rows and 217 columns. One contains the latitudes in each grid cell, the second contains the longitudes, and the third one the values of my represented variable. I define a user variable (type 'i4' in netCDF) called "crs", which contains the EPSG code of interest (EPSG:4326), and the proj4 parameter. The file is successfully written in my disk, but when I open it in QGIS it says that the "CRS is undefined" and then I see the map as blank.

This is how I created the metadata for the netcdf file. Note that the variable "proj_wgs" contains the following string: "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

ds = Dataset(basedir_ou.format(ufile_out), 'w', format='NETCDF4_CLASSIC')

dim_lat = ds.createDimension("x", size=234)
dim_lon = ds.createDimension("y", size=217)
dim_tim = ds.createDimension("time", size=1)
dim_ban = ds.createDimension("bnds", size=2)

ds.createVariable('crs', 'i4')
ds.variables["crs"].grid_mapping_name = 'crs'
ds.variables["crs"].epsg_code = 'EPSG:4326'
ds.variables["crs"].proj4_params = proj_wgs
ds.variables["crs"].spatial_ref = proj_wgs
ds.variables['crs'].valid_min = -200
ds.variables['crs'].valid_max = 200    

ds.createVariable("x", np.float32, ("x",))
ds.variables["x"].grid_mapping = 'crs'
ds.variables["x"].grid_mapping_name = 'latitude_longitude'
ds.variables['x'].valid_min = -200
ds.variables['x'].valid_max = 200

ds.createVariable("y", np.float32, ("y",))
ds.variables["y"].grid_mapping = 'crs'
ds.variables["y"].grid_mapping_name = 'latitude_longitude'
ds.variables['y'].valid_min = -200
ds.variables['y'].valid_max = 200

ds.createVariable("lon", np.float32, ("y", "x"))
ds.variables["lon"].grid_mapping = 'crs'
ds.variables["lon"].grid_mapping_name = 'latitude_longitude'
ds.variables["lon"].units = "degree_east"
ds.variables['lon'].valid_min = -200
ds.variables['lon'].valid_max = 200
ds.variables["lon"][:,:] = proj_lons

ds.createVariable("lat", np.float32, ("y", "x"))
ds.variables["lat"].grid_mapping = 'crs'
ds.variables["lat"].grid_mapping_name = 'latitude_longitude'
ds.variables["lat"].units = 'degree_north'
ds.variables['lat'].valid_min = -200
ds.variables['lat'].valid_max = 200
ds.variables["lat"][:,:] = proj_lats

ds.createVariable("vble", np.float32, ("y", "x"))
ds.variables["vble"].grid_mapping = 'crs'
ds.variables["vble"].grid_mapping_name = 'latitude_longitude'
ds.variables['vble'].valid_min = -100
ds.variables['vble'].valid_max = 100
ds.variables["vble"][:,:] = array_variable

# These variables are commented for simplicity
# ds.createVariable("time", np.float64, ("time",))
# ds.createVariable("date", np.int32, ("time",))
# ds.createVariable("hms", np.int32, ("time"))
# ds.createVariable("time_bnds", np.float64, ("time", "bnds"))
# ds.createVariable("height", np.float32)

By doing "gdalinfo destination_file.nc", I see the following data:

Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Driver: netCDF/Network Common Data Format
Files: none associated
Size is 512, 512
Coordinate System is `'
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:".../output_folder/destination_file.nc":lon
  SUBDATASET_1_DESC=[217x234] lon (32-bit floating-point)
  SUBDATASET_2_NAME=NETCDF:".../output_folder/destination_file.nc":lat
  SUBDATASET_2_DESC=[217x234] lat (32-bit floating-point)
  SUBDATASET_3_NAME=NETCDF:"/output_folder/destination_file.nc":vble
  SUBDATASET_3_DESC=[217x234] vble (32-bit floating-point)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

Thus, it seems that the CRS is not recognized. Then, the size is incorrect, despite I specified the number of rows and columns (234x217), and when I open the metadata in QGIS, it says that the origin of coordinates is "9.96921e+36,9.96921e+36", which is incorrect as well.

Is there a way of correcting for these issues in python-netcdf4?

How can I make QGIS or GDAL to recognize the CRS?

Answer 1

I had the same question and I found this question to be very helpful. Hence, to solve your problem you need to update your code, following the example I include below. Start by studying the crs, lat and lon variables.

res_in_lat, res_in_lon = 0.05, 0.05    # change accordingly
bounds = {"north_lat" : 50,            # change this dic accordingly
          "south_lat" : 40,
          "east_lon" : 50,
          "west_lon" : 40,
          }

with Dataset(basedir_ou.format(ufile_out), 'w', format='NETCDF4_CLASSIC') as ds:

    dim_lat = ds.createDimension("x", size=234)
    dim_lon = ds.createDimension("y", size=217)
    dim_tim = ds.createDimension("time", size=1)
    dim_ban = ds.createDimension("bnds", size=2)

    crs = ds.createVariable('WGS84', 'c')
    crs.spatial_ref = """GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]"""

    lat = ds.createVariable('lat', np.float32, ('dim_lat',))
    lat.units = 'degrees_north'
    lat[:] = np.arange(bounds["north_lat"],bounds["south_lat"],-1*res_in_lat)

    lon = ds.createVariable('lon', np.float32, ('dim_lon',))
    lon.units = 'degrees_east'
    lon[:] = np.arange(bounds["west_lon"],bounds["east_lon"],0.res_in_lon) 

    ds.createVariable("vble", np.float32, ("lat", "lon"))
    vble.grid_mapping = 'WGS84' # the crs variable name
    vble.grid_mapping_name = 'latitude_longitude'
    vble.valid_min = -100
    vble.valid_max = 100
    vble[:] = array_variable