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?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
RRSC NGP
  • 658
  • 1
  • 8
  • 19

4 Answers4

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.

snowman2
  • 7,321
  • 12
  • 29
  • 54
  • 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 04: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 '20 at 09: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 '20 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 '20 at 13:39
  • Nice answer but the library is buggy. It has syntax and import errors! – wondim Jun 15 '21 at 04:30
  • Are you using python 3.7+? It doesn't support older versions of python. – snowman2 Jun 15 '21 at 13:38
  • Hi @snowman2. Can you tell why this code is not working with this file "https://drive.google.com/file/d/1zw8kSPZPRIu_YDxixIlOByPCY19GNMbz/view?usp=sharing" – RRSC NGP May 14 '22 at 10:35
1

You can try replacing the variable et by

import numpy as np

et = np.rot90(et,2)
et = np.flip(et,1)
Rafael sampaio
  • 11
  • 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.

David Bernat
  • 111
  • 2
  • 1
    Mind editing your answer to say that you can use rioxarray.open_rasterio to use rasterio to transpose/rotate the data? – snowman2 Apr 10 '20 at 01: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)
Zoriana
  • 21
  • 4