Set projection for NetCDF4 in Python

Question

I am using the NetCDF4 library for Python but cannot find how to set the projection to WGS84.

import netCDF4 as nc
dsout = nc.Dataset(outfile, 'r')
time = dsout.createDimension('time', 0)
lat = dsout.createDimension('lat', rows)
lat = dsout.createVariable('lat', 'f4', ('lat',))
lat.standard_name = 'latitude'
lat.units = 'degrees_north'
lat.axis = "Y"
lat[:] = lats
lon = dsout.createDimension('lon', cols)
lon = dsout.createVariable('lon', 'f4', ('lon',))
lon.standard_name = 'longitude'
lon.units = 'degrees_east'
lon.axis = "X"
lon[:] = lons
times = dsout.createVariable('time', 'f4', ('time',))
times.standard_name = 'time'
times.long_name = 'time'
times.units = 'hours since 1970年01月01日 00:00:00'
times.calendar = 'gregorian'
acc_precip = dsout.createVariable(
 'acc_precipitation_amount',
 'f4',
 ('time', 'lat', 'lon'),
 zlib=True,
 complevel=9,
 least_significant_digit=1,
 fill_value=-9999
)
acc_precip.standard_name = 'acc_precipitation_amount'
acc_precip.units = 'mm'

• Gurminder Bharani: Do you have access to ArcMap or GDAL?
• Jens de Bruijn: GDAL, yes. However, I would prefer a solution within the same library :)
• Gurminder Bharani: Ok, are you trying to export this raster data into tiff or img with WGS84 projection?
• Jens de Bruijn: No, I am trying to add the WGS84 projection to the NetCDF files. In other NetCDF files that I did not create it says it has the WGS projection. When I create the files myself using the above method, I do not know how to add it. Hence, my question.

Answer 1

I don't know much at all about writing NetCDF files, but the following allowed me to write a NetCDF file that GDAL recognised the CRS.

import numpy as np
import netCDF4 as nc
outfile = r'test.nc'
lats = tuple(range(-90, 90))
lons = tuple(range(0, 180))
cols = len(lons)
rows = len(lats)
dsout = nc.Dataset(outfile, 'w', clobber=True)
time = dsout.createDimension('time', 0)
lat = dsout.createDimension('lat', rows)
lat = dsout.createVariable('lat', 'f4', ('lat',))
lat.standard_name = 'latitude'
lat.units = 'degrees_north'
lat.axis = "Y"
lat[:] = lats
lon = dsout.createDimension('lon', cols)
lon = dsout.createVariable('lon', 'f4', ('lon',))
lon.standard_name = 'longitude'
lon.units = 'degrees_east'
lon.axis = "X"
lon[:] = lons
times = dsout.createVariable('time', 'f4', ('time',))
times.standard_name = 'time'
times.long_name = 'time'
times.units = 'hours since 1970年01月01日 00:00:00'
times.calendar = 'gregorian'
acc_precip = dsout.createVariable(
 'acc_precipitation_amount',
 'f4',
 ('time', 'lat', 'lon'),
 zlib=True,
 complevel=9,
 least_significant_digit=1,
 fill_value=-9999
)
acc_precip[:] = np.ones((1,180,180))
acc_precip.standard_name = 'acc_precipitation_amount'
acc_precip.units = 'mm'
acc_precip.setncattr('grid_mapping', 'spatial_ref')
crs = dsout.createVariable('spatial_ref', 'i4')
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.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]'

gdalinfo output:

$ gdalinfo netcdf:test.nc
Warning 1: NetCDF driver detected file type=5, but libnetcdf detected type=3
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Driver: netCDF/Network Common Data Format
Files: none associated
Size is 180, 180
Coordinate System is:
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.0174532925199433,
 AUTHORITY["EPSG","9122"]],
 AUTHORITY["EPSG","4326"]]
Origin = (-0.500000000000000,89.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
 acc_precipitation_amount#grid_mapping=spatial_ref
 acc_precipitation_amount#least_significant_digit=1
 acc_precipitation_amount#standard_name=acc_precipitation_amount
 acc_precipitation_amount#units=mm
 acc_precipitation_amount#_FillValue=-9999
 lat#axis=Y
 lat#standard_name=latitude
 lat#units=degrees_north
 lon#axis=X
 lon#standard_name=longitude
 lon#units=degrees_east
 NETCDF_DIM_EXTRA={time}
 NETCDF_DIM_time_DEF={1,5}
 NETCDF_DIM_time_VALUES=9.96921e+36
 spatial_ref#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.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]
 time#calendar=gregorian
 time#long_name=time
 time#standard_name=time
 time#units=hours since 1970年01月01日 00:00:00
Corner Coordinates:
Upper Left ( -0.5000000, 89.5000000) ( 0d30' 0.00"W, 89d30' 0.00"N)
Lower Left ( -0.5000000, -90.5000000) ( 0d30' 0.00"W, 90d30' 0.00"S)
Upper Right ( 179.5000000, 89.5000000) (179d30' 0.00"E, 89d30' 0.00"N)
Lower Right ( 179.5000000, -90.5000000) (179d30' 0.00"E, 90d30' 0.00"S)
Center ( 89.5000000, -0.5000000) ( 89d30' 0.00"E, 0d30' 0.00"S)
Band 1 Block=180x1 Type=Float32, ColorInterp=Undefined
 NoData Value=-9999
 Unit Type: mm
 Metadata:
 grid_mapping=spatial_ref
 least_significant_digit=1
 NETCDF_DIM_time=9.96921e+36
 NETCDF_VARNAME=acc_precipitation_amount
 standard_name=acc_precipitation_amount
 units=mm
 _FillValue=-9999