How to convert a netCDF4 file to a geoTiff

Question

I'm currently trying to get Tropomi data in geoTiff format. I downloaded some data in netCDF4 format. This way I obtain three numpy arrays. one with latitude coordinates, one with longitude coordinates and one with carbon-mono-oxide values.

So I have a matrix with values for my raster and of each value I know the longitude and latitude of that respective value.

With this information how can I construct a georeferenced raster?

I read in the data as follows import netCDF4 from netCDF4 import Dataset import numpy as np

file = '/home/daniel/Downloads/S5P_NRTI_L2__CO_____20190430T171319_20190430T171819_08006_01_010301_20190430T175151.nc'

rootgrp = Dataset(file, "r",format="NETCDF4")

lat = rootgrp.groups['PRODUCT']['latitude'][:] 
lon = rootgrp.groups['PRODUCT']['longitude'][:]
carbon = rootgrp.groups['PRODUCT']['carbonmonoxide_total_column'][:]

obtaining 3 matrices with shape (1,290,215)

Now I would like to convert this to a Mercator projected geoTIFF, but I do not know how to go about it.

Possible duplicate of [Convert NetCDF (.nc) to GEOTIFF](https://stackoverflow.com/questions/52046282/convert-netcdf-nc-to-geotiff) — ClimateUnboxed, May 01 '19 at 15:03

score 2 · Answer 1 · answered May 02 '19 at 12:03

the gdal_translate option seems to work. But here is an alternative explicit way I did it.

#importing packages
import numpy as np
from scipy import interpolate
from netCDF4 import Dataset
from shapely.geometry import Point
import geopandas as gpd
from geopy.distance import geodesic
import rasterio
import matplotlib.pyplot as plt

#load data 
file = '/home/daniel/Ellipsis/db/downloaded/rawtropomi/S5P_NRTI_L2__CO_____20190430T171319_20190430T171819_08006_01_010301_20190430T175151.nc'
rootgrp = Dataset(file, "r",format="NETCDF4")
lat = rootgrp.groups['PRODUCT']['latitude'][:]
lon = rootgrp.groups['PRODUCT']['longitude'][:]
carbon = rootgrp.groups['PRODUCT']['carbonmonoxide_total_column'][:]
carbon = carbon.filled(0)
lat = lat.filled(-1000)
lon = lon.filled(-1000)

carbon = carbon.flatten()
lat = lat.flatten()
lon = lon.flatten()

#calculate the real distance between corners and get the widht and height in pixels assuming you want a pixel resolution of at least 7 by 7 kilometers
w = max(geodesic((min(lat),max(lon)), (min(lat),min(lon))).meters/7000 , geodesic((max(lat),max(lon)), (max(lat),min(lon))).meters/14000)
h = geodesic((min(lat),max(lon)), (max(lat),max(lon))).meters/14000

# create a geopandas with as its rows the latitude, longitude an the measrument values. transfrom it to the webmercator projection (or projection of your choosing)
points = [Point(xy) for xy in zip(lon, lat)]
crs = {'init': 'epsg:4326'}
data = gpd.GeoDataFrame({'value':carbon}, crs=crs, geometry=points)
data = data.to_crs({'init': 'epsg:3395'})
data['lon'] = data.bounds['maxx'].values
data['lat'] = data.bounds['maxy'].values

#make grid of coordinates. You nee de calculate the coordinate of each pixel in the desired raster
minlon = min(data['lon'])
maxlon = max(data['lon'])
minlat = min(data['lat'])
maxlat = max(data['lat'])

lon_list = np.arange(minlon, maxlon, (maxlon-minlon)/w )
lat_list = np.arange(minlat, maxlat, (maxlat-minlat)/h)

lon_2d, lat_2d = np.meshgrid(lon_list, lat_list)



#use the values in the geopandas dataframe to interpolate values int the coordinate raster
r = interpolate.griddata(points = (data['lon'].values,data['lat'].values), values = data['value'].values, xi = (lon_2d, lat_2d))
r = np.flip(r, axis = 0)

#check result
plt.imshow(r)


#save raster
transform = rasterio.transform.from_bounds(south = minlat, east = maxlon, north =     maxlat, west = minlon, width = r.shape[1], height = r.shape[2]   )
file_out = 'test.tiff'
new_dataset = rasterio.open(file_out , 'w', driver='Gtiff', compress='lzw',
                                    height = r.shape[1], width = r.shape[2],
                                    count= r.shape[0], dtype=str( r.dtype),
                                    crs=   data.crs,
                                    transform= transform)
new_dataset.write(r)
new_dataset.close()

Your method seems useful. But I can't work out why you flatten your latitude and longitude from the netcdf: shouldn't they be 1D arrays anyway? — Pauli, Nov 20 '20 at 08:58
unfortunately whgen attempting this, my `r` only has two dimensions (I change `height` to `r.shape[0]` and width to `r,shape[1]` and count = 1 in `rasterio.open`) as it returns a tuple dimension error. However i now have the error `Source shape is inconsistent with given indexes`. Should `r` have 3 dimensions? — Sam, Sep 13 '21 at 10:53

score 1 · Answer 2 · answered May 01 '19 at 15:02

1

I would suggest looking at this answer here using gdal_translate:

Convert NetCDF (.nc) to GEOTIFF

gdal_translate -of GTiff file.nc test.tiff

answered May 01 '19 at 15:02

ClimateUnboxed

7,106
3
41
86

How to convert a netCDF4 file to a geoTiff

2 Answers2