Converting NetCDF dataset array to GeoTiff using rasterio Python

Question

I have read a NetCDF file using the netCDF4 library and then read one of its datasets ("Evapotranspiration") into a variable (variable contains array) using the following code. Subsequently now I am trying to convert this array into a GeoTIFF file using rasterio. However, the resulting GeoTIFF is appearing to be rotated by 90 Degrees when I am opening it in QGIS. Following is my code:

from netCDF4 import Dataset
import rasterio
from rasterio.transform import from_origin
nc = Dataset(r'D:\Weather data\et_01012018.nc','r')
 
lat = nc.variables['latitude'][:]
lon = nc.variables['longitude'][:]
et = nc.variables['Evapotranspiration'][:]
transform = from_origin(68.175 , 37.025 , 0.05, 0.05)
profile = {'driver': 'GTiff', 'height': et.shape[1], 'width': et.shape[0], 'count': 1, 'dtype': str(et.dtype), 'transform': transform}
with rasterio.open(r'D:\Weather data\test.tif', 'w', crs='EPSG:4326', **profile) as dst:
 dst.write(et,1)

Further I also tried GDAL to implement the same but no success as of now. I am getting same results i.e. the raster file is rotated by 90 degrees in clockwise direction. Following is the code implemented by using GDAL:

import gdal, osr
from netCDF4 import Dataset
nc = Dataset(r'D:\Weather data\et_01012018.nc','r')
lat = nc.variables['latitude'][:]
lon = nc.variables['longitude'][:]
et = nc.variables['Evapotranspiration'][:]
# reverse array so the tif looks like the array
et_T = et[::-1] 
cols = et_T.shape[1]
rows = et_T.shape[0]
# Origin should be in Longitude-Latitude form
originX = 79
originY = 21
driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create(r'D:\Weather data\test.tif', cols, rows, 1, gdal.GDT_Float64)
outRaster.SetGeoTransform((originX, 0.05, 0, originY, 0, 0.05))
outband = outRaster.GetRasterBand(1)
outband.WriteArray(et_T)
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromEPSG(4326)
outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()
outRaster = None

Can you help me to resolve the issue?

• bugmenot123: Does this also happen if you use a GeoTIFF as input data?
• RRSC NGP: @Bugmenot123 I haven't tried it....will try and see....

Answer 1

I would recommend looking into rioxarray for your dataset.

You can open your dataset like so:

import rioxarray
import xarray
xds = xarray.open_dataset('D:\Weather data\et_01012018.nc')

If your CRS is not discovered, you should be able to add it like so:

xds.rio.write_crs("epsg:4326", inplace=True)

Then, you should be able to create a geotiff from the Evapotranspiration like so:

xds["Evapotranspiration"].rio.to_raster('D:\Weather data\test.tif')

If this does not produce the correct results, I would be interested in learning more about your input file data.

• RRSC NGP: Hi @snowman2....thanks for your reply. Ur solution is working fine for me. In some of the NC files "InvalidDimensionOrder" error is coming which I think is the issue while creating the NC files. As a solution I am transposing the dataset and the issue is getting resolved.
• snowman2: Great, I am glad to hear that this solution works out of the box.
• MCMZL: @snowman2 it looks to me that it will only work if your longitude and latitude coordinates are already defined with an affine transform, am I correct? I am looking for an example with irregular grid
• snowman2: This example only works with a regular grid. If you have an irregular grid, that usually means that your data was projected into another projection. The solution in that case is to convert your data back into the original projection with a regular grid.
• snowman2: If that is not the case, then you will need to resample your data onto a regular grid for this to work.

Answer 2

You can try replacing the variable et by

import numpy as np
et = np.rot90(et,2)
et = np.flip(et,1)

Answer 3

For anyone arriving here from Google, it does appear that rasterio does read netCDF files such that the data is transposed and rotated from the typical reading of GeoTiffs.

• snowman2: Mind editing your answer to say that you can use rioxarray.open_rasterio to use rasterio to transpose/rotate the data?

Answer 4

You can read nc file directly with rasterio.

import rasterio as rio
var = "Evapotranspiration"
reader = rio.open(r"netcdf:Path\to\nc_file\nc_file.nc:"+var)
prof = reader.profile
prof.update(driver='GTiff', crs="EPSG:4326",)
array = reader.read(1)
with rio.open(rf"Path\to\new_tiff_file\new_tif.tif", "w", **prof) as dst:
 dst.write(array,1)