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
Best Answer
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
try this:
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.