Determine the geotransformation to convert a netCDF to geotiff

Question

I have a netCDF file I need to convert to a geotiff. There are three subdatasets in the netCDF file: snow, lat and lon. I am accessing these files like this:

import gdal
import netCDF4 as nc

ncfile = nc.Dataset("C:/path_to_file/file.nc", "r")

snow = ncfile.variables["snc"]
lat = ncfile.variables["lat"]
lon = ncfile.variables["lon"]

The metadata says the latitude and longitude are in degrees, so I want to project to WGS.

The minimum and maximum latitude and longitudes found by:

xmin, ymin, xmax, ymax = [lon.min(), lat.min(), lon.max(), lat.max()]

returns:

(0.0, -87.863801343710804, 357.1875, 87.863801343710804)

since my longitude goes from 0 to 360, and it should go from -180 to 180 I am adjusting it like so:

lon = ((lon+180) % 360) - 180

I now determine my geotransform:

nx = len(lon)
ny = len(lat)
xmin, ymin, xmax, ymax = [lon.min(), lat.min(), lon.max(), lat.max()]
xres = (xmax - xmin) / float(nx)
yres = (ymax - ymin) / float(ny)
geotransform = (xmin, xres, 0, ymax, 0, -yres)

and try to create my tif:

dst_ds = gdal.GetDriverByName('GTiff').Create('output.tif', ny, nx, 1, gdal.GDT_Float32)


dst_ds.SetGeoTransform(geotransform)    # specify coords
srs = osr.SpatialReference()            # establish encoding
srs.ImportFromEPSG(3857)                # WGS84 lat/long
dst_ds.SetProjection(srs.ExportToWkt()) # export coords to file
dst_ds.GetRasterBand(1).WriteArray(snow)   # write r-band to the raster
dst_ds.FlushCache()                     # write to disk
dst_ds = None     

I followed this unanswered question to get this far, https://stackoverflow.com/questions/43254024/extracting-specific-netcdf-info-and-converting-to-geotiff-in-python, but my raster won't even output.

I instead get this error:

Traceback (most recent call last):

  File "<ipython-input-10-0384d872ba66>", line 15, in <module>
    dst_ds.GetRasterBand(1).WriteArray(snow)   # write r-band to the raster

  File "C:\Users\spotter\AppData\Local\Continuum\Anaconda_64\lib\site-packages\osgeo\gdal.py", line 2332, in WriteArray
    callback_data = callback_data )

  File "C:\Users\spotter\AppData\Local\Continuum\Anaconda_64\lib\site-packages\osgeo\gdal_array.py", line 365, in BandWriteArray
    raise ValueError("array larger than output file, or offset off edge")

ValueError: array larger than output file, or offset off edge

Answer 1

For me, it looks like the "lat"/"lon" bands do not have the same number of pixels as the "snow" band. I experienced a similar behaviour with Sentinel-3 data. Instead of

nx = len(lon)
ny = len(lat)

try this:

nx = ncfile.dimensions['columns'].size
ny = ncfile.dimensions['rows'].size

Maybe your netCDF has different names for the dimensions, you should look into the metadata for this or just try it out in the Python console.

Additionally, EPSG:3857 is not WGS84 lat/lon, but WGS84 Web Mercator. You are looking for EPSG:4326. However, when using Sentinel-3 data, I'm still getting an inacceptable spatial offset, but this might not be the case with your data.

Answer 2

This code requires minor changes and it will work. Instead of:

dst_ds = gdal.GetDriverByName('GTiff').Create('output.tif', ny, nx, 1, gdal.GDT_Float32)

Use:

dst_ds = gdal.GetDriverByName('GTiff').Create('output.tif', nx, ny, 1, gdal.GDT_Float32)

Answer 3

I also want to comment on how the variables are being read from the netcdf file. read varaiblaes as follows:

import gdal
import netCDF4 as nc

ncfile = nc.Dataset("C:/path_to_file/file.nc", "r")

snow = ncfile.variables["snc"][:,:]
lat = ncfile.variables["lat"][:]
lon = ncfile.variables["lon"][:]