Using python to write irregularly spaced, georefernced 2D array to netcdf

user_13204380 0
Viewed 269

I have an array (1500 x 2500) of data that I want to write to a netcdf file, but the lat/lon grid is not regularly spaced. I can create a netcdf file with data in it, but the file doesn't appear to be readable and/or in the proper format. For example, when I try to read it with the NOAA Weather and Climate Toolkit, I get an error that states the file was scanned as a 'gridded' file, but no grids were found. Here is my code:

from netCDF4 import Dataset

ndvi_nc = Dataset("ndvi_out.nc", "w", format="NETCDF4")

print(ndvi_nc.data_model)

time = ndvi_nc.createDimension("time", 1)
lat = ndvi_nc.createDimension("lat", 3750000)
lon = ndvi_nc.createDimension("lon", 3750000)
data = ndvi_nc.createDimension("data", 3750000)

print(lats.shape)
print(lons.shape)
print(NDVInew.shape)

times = ndvi_nc.createVariable("time","f8",("time",))
latitudes = ndvi_nc.createVariable("lat","f4",("lat",))
longitudes = ndvi_nc.createVariable("lon","f4",("lon",))
ndvi = ndvi_nc.createVariable("ndvi","f4",("data",))

ndvi_nc.description = "NDVI dataset"
ndvi_nc.source = "netCDF4 python module tutorial"
latitudes.units = "degrees north"
longitudes.units = "degrees west"
ndvi.units = "dimensionless"
times.units = "seconds since 1970-1-1"

flatlats = lats.flatten()
flatlons = lons.flatten()
flatndvi = NDVInew.flatten()
print(flatndvi.shape)

latitudes[:] = flatlats
longitudes[:] = flatlons
ndvi[:] = flatndvi
times = 1303939211

print(ndvi)
print(latitudes)
print(longitudes)

ndvi_nc.close()

And here is the output:

NETCDF4
(1500, 2500)
(1500, 2500)
(1500, 2500)
(3750000,)
<class 'netCDF4._netCDF4.Variable'>
float32 ndvi(data)
    units: dimensionless
unlimited dimensions: 
current shape = (3750000,)
filling on, default _FillValue of 9.969209968386869e+36 used
<class 'netCDF4._netCDF4.Variable'>
float32 lat(lat)
    units: degrees north
unlimited dimensions: 
current shape = (3750000,)
filling on, default _FillValue of 9.969209968386869e+36 used
<class 'netCDF4._netCDF4.Variable'>
float32 lon(lon)
    units: degrees west
unlimited dimensions: 
current shape = (3750000,)
filling on, default _FillValue of 9.969209968386869e+36 used

Is there a better way to be writing the data to netcdf?

Thanks,

-- Matt

netcdf output grid python
1 Answers

You don't need to flatten everything, it's straightforward to write a grid to netCDF with uneven spacing. You create two dimensions x/y, and then lon, lat, and ndvi are all 2D fields on those dimensions. So something like this worked for me:

from netCDF4 import Dataset
ndvi_nc = Dataset("ndvi.nc", "w")

lats = np.random.rand(1500, 2500)
lons = np.random.rand(1500, 2500)
NDVInew = np.random.rand(1500, 2500)

time = ndvi_nc.createDimension("time", 1)
lat = ndvi_nc.createDimension("y", 1500)
lon = ndvi_nc.createDimension("x", 2500)

times = ndvi_nc.createVariable("time", "f8", ("time",))
latitudes = ndvi_nc.createVariable("lat", "f4", ("y", "x"))
longitudes = ndvi_nc.createVariable("lon", "f4", ("y", "x"))
ndvi = ndvi_nc.createVariable("ndvi", "f4", ("y", "x"))

ndvi_nc.description = "NDVI dataset"
ndvi_nc.source = "netCDF4 python module tutorial"

latitudes.units = "degrees north"
longitudes.units = "degrees west"
ndvi.units = "dimensionless"
times.units = "seconds since 1970-1-1"

flatlats = lats.flatten()
flatlons = lons.flatten()
flatndvi = NDVInew.flatten()

latitudes[:] = lats
longitudes[:] = lons
ndvi[:] = NDVInew
times[:] = 1303939211

ndvi_nc.close()
user_1193140
Related