1. This site uses cookies. By continuing to use this site, you are agreeing to our use of cookies. Learn More.

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

Discussion in 'Geography' started by Irene, Oct 8, 2018.

  1. Irene

    Irene Guest

    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?

    Login To add answer/comment
     

Share This Page