stcorp / harp

Data harmonization toolset for scientific earth observation data
http://stcorp.github.io/harp/doc/html/index.html
BSD 3-Clause "New" or "Revised" License
55 stars 19 forks source link

TROPOMI-Sentinel-5P_L2-TROPOMAER #300

Closed NoOne-play closed 1 year ago

NoOne-play commented 1 year ago

I want to merge this file but getting unsupported product error. How to address this issue. i try other methods like cdo, xarray, none of them works, I think this is because this has it has two groups GEODATA and SCIDATA. how do i spatially merge multiple to this nc file. Any insight please.

svniemeijer commented 1 year ago

In HARP we currently only support the S5P products that are provided by the Copernicus Data Space Environment and by S5P-PAL. If I understand correctly, the TROPOMAER product is a NASA-specific product. I can check if we can add support for this in HARP itself. It is this product, right?

Note that you can also always create a HARP converter yourself, by e.g. reading the data in python using the netcdf library, constructing a harp.Product() with the variables that you need (see the examples for how this can be done), and exporting this to a HARP compatible file with harp.export_product.

NoOne-play commented 1 year ago

Thank you for your reply. Yes I am using NASA-specific product that you mention. I tried creating HARP converter as you suggested with the help of this. This is my code, kindly take a look at it it is giving me TypeError: iteration over a 0-d array. particularly due to this line times = np.array([(t - time_origin).total_seconds() / (24 60 60) for t in times]). Thanks again for your help.

This is my code import netCDF4 import numpy as np from datetime import datetime, timedelta import harp

nc_file = "C:/Users/shekh/OneDrive/Desktop/New folder/test_j/TROPOMI-Sentinel-5P_L2-TROPOMAER_2019m0508t071233-o08114_v01-2021m1011t031834.nc" dataset = netCDF4.Dataset(nc_file)

times = dataset['GEODATA']['time'][:] latitude_bounds = dataset['GEODATA']['latitude_bounds'][:] longitude_bounds = dataset['GEODATA']['longitude_bounds'][:] FinalAerosolOpticalDepth = dataset['SCIDATA']['FinalAerosolOpticalDepth'][:]

product = harp.Product()

time_origin = datetime(2010, 1, 1) times = np.array([(t - time_origin).total_seconds() / (24 60 60) for t in times]) product.datetime = harp.Variable(times, ["time"]) product.datetime.unit = "days since 2010-01-01" # Adjust the format to match your time origin

product.latitude_bounds = harp.Variable(latitude_bounds, ["time", None]) product.latitude_bounds.unit = "degree_north"

product.longitude_bounds = harp.Variable(longitude_bounds, ["time", None]) product.longitude_bounds.unit = "degree_east"

product.FinalAerosolOpticalDepth = harp.Variable(FinalAerosolOpticalDepth, ["time"]) product.FinalAerosolOpticalDepth.unit = "unit"

harp.export_product(product, "myproduct.nc")

svniemeijer commented 1 year ago

These questions are actually better asked on the forum.

I've updated the script to provide a proper HARP export, including only valid pixels. This should get you going.

import netCDF4
import numpy as np
from datetime import datetime, timedelta
import harp

nc_file = "TROPOMI-Sentinel-5P_L2-TROPOMAER_2019m0508t071233-o08114_v01-2021m1011t031834.nc"
dataset = netCDF4.Dataset(nc_file)

times = dataset['GEODATA']['delta_time'][:]
latitude_bounds = dataset['GEODATA']['latitude_bounds'][:]
longitude_bounds = dataset['GEODATA']['longitude_bounds'][:]
FinalAerosolOpticalDepth = dataset['SCIDATA']['FinalAerosolOpticalDepth'][:]
shape = dataset['SCIDATA']['FinalAerosolOpticalDepth'].shape

product = harp.Product()

times = harp.convert_unit(dataset['GEODATA']['delta_time'].units, "seconds since 2010-01-01 00:00:00", times)
times = np.repeat(times, shape[1])
product.datetime = harp.Variable(times, ["time"])
product.datetime.unit = "seconds since 2010-01-01"

latitude_bounds = np.reshape(latitude_bounds, (shape[0] * shape[1], 4))
product.latitude_bounds = harp.Variable(latitude_bounds, ["time", None])
product.latitude_bounds.unit = "degree_north"

longitude_bounds = np.reshape(longitude_bounds, (shape[0] * shape[1], 4))
product.longitude_bounds = harp.Variable(longitude_bounds, ["time", None])
product.longitude_bounds.unit = "degree_east"

FinalAerosolOpticalDepth = np.reshape(FinalAerosolOpticalDepth, (shape[0] * shape[1], shape[2]))
# pick a specific wavelength
FinalAerosolOpticalDepth = FinalAerosolOpticalDepth[:,0]
product.aerosol_optical_depth = harp.Variable(FinalAerosolOpticalDepth, ["time"])
product.aerosol_optical_depth.unit = "1"
product.aerosol_optical_depth.valid_min = 0
product.aerosol_optical_depth.valid_max = 4

harp.export_product(product, "myproduct.nc", operations="valid(aerosol_optical_depth)")
NoOne-play commented 1 year ago

Thank you very much.

NoOne-play commented 1 year ago

Solved,