5

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

Kindly help me in resolving the issue.

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

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.set_crs("epsg:4326")

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.

| improve this answer | |
  • 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 Jul 22 '19 at 4:26
  • Great, I am glad to hear that this solution works out of the box. – snowman2 Jul 23 '19 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 Aug 17 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 Aug 17 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 Aug 17 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)
| improve this answer | |
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.

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

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.