12

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?

2
  • Does this also happen if you use a GeoTIFF as input data? Commented May 20, 2019 at 13:18
  • @Bugmenot123 I haven't tried it....will try and see....
    – RRSC NGP
    Commented May 20, 2019 at 13:44

4 Answers 4

20

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.

8
  • 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.
    – RRSC NGP
    Commented Jul 22, 2019 at 4:26
  • Great, I am glad to hear that this solution works out of the box.
    – snowman2
    Commented Jul 23, 2019 at 12:12
  • @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
    – MCMZL
    Commented Aug 17, 2020 at 9:43
  • 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
    Commented Aug 17, 2020 at 13:39
  • If that is not the case, then you will need to resample your data onto a regular grid for this to work.
    – snowman2
    Commented Aug 17, 2020 at 13:39
1

You can try replacing the variable et by

import numpy as np

et = np.rot90(et,2)
et = np.flip(et,1)
1

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.

1
  • 1
    Mind editing your answer to say that you can use rioxarray.open_rasterio to use rasterio to transpose/rotate the data?
    – snowman2
    Commented Apr 10, 2020 at 1:21
0

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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