Issue with converting Sentinel-5P NetCDF file to GeoTIFF file

[rrsc1234]

I am trying to convert Sentinel-5P NetCDF file (available at https://drive.google.com/file/d/1dJdhQspdI3p5YyyiBBiNxUPCCOU91QrU/view?usp=sharing) to GeoTIFF using the following code:

import xarray
import rioxarray

netcdf_fname = r'C:\Users\User\Desktop\NC Files\CH4_1.nc'

xds = xarray.open_dataset(netcdf_fname, group='PRODUCT')

param = xds.methane_mixing_ratio
param.rio.set_spatial_dims(x_dim="longitude", y_dim="latitude", inplace=True)
param.rio.write_crs("EPSG:4326", inplace=True)

param.rio.to_raster(r'C:\Users\User\Desktop\NC Files\ch4_raster.tif')

After running the above code I am getting the following error: MissingSpatialDimensionError: x dimension (longitude) not found. Data variable: methane_mixing_ratio

How can I resolve this issue?

[snowman2]

How can I resolve this issue?

I recommend selecting the subset of data variables that have data on a geospatial grid before exporting. The error is coming from an attempt to export a non-gridded variable to a raster.

[rrsc1234]

How can I resolve this issue?

I recommend selecting the subset of data variables that have data on a geospatial grid before exporting. The error is coming from an attempt to export a non-gridded variable to a raster.

I tried with other parameters available in the dataset and getting the same result as mentioned earlier:

import xarray
import rioxarray

netcdf_fname = r'C:\Users\User\Desktop\NC Files\CH4_1.nc'

xds = xarray.open_dataset(netcdf_fname, group='PRODUCT')

nc_vars = list(xds._variables.keys())
# scanline
# ground_pixel
# time
# corner
# layer
# level
# delta_time
# time_utc
# qa_value
# latitude
# longitude
# methane_mixing_ratio
# methane_mixing_ratio_precision
# methane_mixing_ratio_bias_corrected

latitude = xds['latitude'][0].values  ## Array of Float32 (4173, 215) 
longitude = xds['longitude'][0].values  ## Array of Float32 (4173, 215) 

methane_mixing_ratio = xds['methane_mixing_ratio'][0].values  ## Array of Float32 (4173, 215)
methane_mixing_ratio_precision = xds['methane_mixing_ratio_precision'][0].values  ## Array of Float32 (4173, 215)
methane_mixing_ratio_bias_corrected = xds['methane_mixing_ratio_bias_corrected'][0].values  ## Array of Float32 (4173, 215)

param = xds.methane_mixing_ratio
# param = xds.methane_mixing_ratio_precision ## Getting the same result as of above parameter
# param = xds.methane_mixing_ratio_bias_corrected ## Getting the same result as of above parameter

param.rio.set_spatial_dims(x_dim="longitude", y_dim="latitude", inplace=True)
param.rio.write_crs("EPSG:4326", inplace=True)

param.rio.to_raster(r'C:\Users\User\Desktop\NC Files\ch4_raster.tif')

Kindly suggest.

[snowman2]

I tried with other parameters available in the dataset and getting the same result as mentioned earlier

It sounds like a valid grid doesn't exist in the dataset.

[rrsc1234]

It sounds like a valid grid doesn't exist in the dataset.

I had tried the following code:

import xarray as xr
# import numpy as np

nc_file_path = r'C:\Users\User\Desktop\NC Files\CH4_1.nc'

ncfile = xr.open_dataset(nc_file_path, group='PRODUCT')

## Extract the variable and converting that into list
param = list(ncfile['methane_mixing_ratio_precision'].values[0].ravel())
# param_min = np.nanmin(param)  ## 1.0205439
# param_max = np.nanmax(param)  ## 36.001057

## Getting the data type of the parameter
data_type = ncfile['methane_mixing_ratio_precision'].values[0].dtype  ## 'float32'

## Extract the latitude and longitud and converting into list
latitude = list(ncfile['latitude'].values[0].ravel())
longitude = list(ncfile['longitude'].values[0].ravel())

# latitude_min = np.nanmin(latitude)  ## -82.8837
# latitude_max = np.nanmax(latitude)  ## 89.968956
# longitude_min = np.nanmin(longitude)  ## -179.83487
# longitude_max = np.nanmax(longitude)  ## 179.93947

## Making a Geodataframe from the list of param, latitude, longitude 
import geopandas as gpd

gdf = gpd.GeoDataFrame(columns = ['Param', 'Lat', 'Lon'])

gdf['Param'] = param
gdf['Lat'] = latitude
gdf['Lon'] = longitude

## Creating a geometry column from latitude and longitude column
gdf['geometry'] = gpd.points_from_xy(gdf.Lon, gdf.Lat, crs="EPSG:4326")
gdf_valid_data = gdf[gdf['Param'].notna()]

# gdf = gdf.to_crs(epsg = 4326)
gdf_valid_data = gdf_valid_data.to_crs(epsg = 4326)

# gdf.plot()
# gdf_valid_data.plot()

gdf_valid_data.to_file(r'C:\Users\User\Desktop\NC Files\methane_mixing_ratio_precision_points.shp')

After running the above code I am able to get point shapefile where the parameter has valid data https://drive.google.com/file/d/1bPE3ssiN4upyStThPuZMHew-Myzngwos/view?usp=sharing.

The above shapefile was verified by visualizing the output of the below code:

import xarray
import matplotlib.pyplot as plt

netcdf_fname = r'C:\Users\User\Desktop\NC Files\CH4_1.nc'

xds = xarray.open_dataset(netcdf_fname, group='PRODUCT')

plt.figure(figsize=(14,8))
ax = plt.axes()

xds.methane_mixing_ratio_precision[0].plot.pcolormesh(ax=ax, x='longitude', 
                                                 y='latitude',
                                                 add_colorbar=True, cmap='jet');

The output of the above code is https://drive.google.com/file/d/1iisKq263bX7wPo8eTWZ9zTyVbUQor-rG/view?usp=sharing.

Now I am stuck with how to convert this point shapefile to polygon shapefile which then I can convert to raster.

[snowman2]

209 might be a helpful reference.

[rrsc1234]

209 might be a helpful reference.

I had tried your suggestion using geocube library:

import geopandas as gpd

gdf = gpd.read_file(r'C:\Users\User\Desktop\NC Files\methane_mixing_ratio_precision_points.shp')

gdf.crs
'''
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
'''

from geocube.api.core import make_geocube
from geocube.rasterize import rasterize_points_griddata

geo_grid = make_geocube(
    vector_data=gdf,
    measurements=['Param'],
    resolution=(-0.03153153153153153, 0.04954954954954955),
    rasterize_function=rasterize_points_griddata,
)

geo_grid.Param.rio.to_raster(r'C:\Users\User\Desktop\NC Files\methane_mixing_ratio_precision_grids.tif')

The GeoTIFF created shows absurd result https://drive.google.com/file/d/1qMKg_ZJOzFAjJRwVNF1aBKzxeLlU7CyW/view?usp=sharing.

Can you please suggest how to resolve this issue.

[snowman2]

This may be helpful: https://github.com/bopen/xarray-sentinel

[rrsc1234]

This may be helpful: https://github.com/bopen/xarray-sentinel

I think this is related to sentinel 1 data processing.

[snowman2]

This may be helpful: https://github.com/bopen/xarray-sentinel

I think this is related to sentinel 1 data processing.

The code could be a good reference as there might be similarities. The developers might be able to offer suggestions. And, if you come up with a solution for your data, you may be able to contribute it there.