4

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
  • Just a guess. You may need to convert to geoTIFF leaving the 0-360 longitude range, then manipulate the geoTIFF to move the >180 portion to the "left". – mkennedy Nov 13 '17 at 20:34
  • without adjustment the same error is returned. – Stefano Potter Nov 13 '17 at 20:54
  • @mkennedy I output a raster correctly using the arcpy tool, but now I still need to adjust to the left like you say. Do you have any ideas how this could be done? – Stefano Potter Nov 13 '17 at 21:28
  • I am finding Q&A to change the range using R and GDAL but not ArcGIS. Here's the GDAL one: gis.stackexchange.com/questions/37790/… – mkennedy Nov 13 '17 at 23:48
2

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.

0

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)
0

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"][:]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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