7

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'
Improve this question
4
  • Do you have access to ArcMap or GDAL?
    – Gurminder Bharani
    Commented Feb 27, 2017 at 15:33
  • GDAL, yes. However, I would prefer a solution within the same library :)
    – Jens de Bruijn
    Commented Feb 27, 2017 at 15:57
  • Ok, are you trying to export this raster data into tiff or img with WGS84 projection?
    – Gurminder Bharani
    Commented Feb 27, 2017 at 18:38
  • 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.
    – Jens de Bruijn
    Commented Feb 27, 2017 at 21:32

1 Answer 1

Reset to default
4

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
Improve this answer
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.