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?

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.

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.

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)