# Determine the geotransformation to convert a netCDF to geotiff

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.

``````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". Commented Nov 13, 2017 at 20:34
• without adjustment the same error is returned. Commented Nov 13, 2017 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? Commented Nov 13, 2017 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/… Commented Nov 13, 2017 at 23:48

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.

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

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