Problems with Harp on Pyhton

HARP
,
1

Hello!

I have some problems with harp.import_product():

  1. The output of the command below is 435x494 but I specify 436x495. Even if I select a wider area when first filtering this is still giving me the same dimensions.

  2. After creating the grid in the command below with bin_spatial(…), I use derive(…) to the get the center latitude and longitude per pixel. However, once I load the object with xarray I only have as coordinates the time. Can you help me to get the expected coordinates?

  3. Is there any way to explore the object after applying harp.import_product() within python other than harp.export_product(…) and then loading it with xarray?

command:
operations:
tropospheric_NO2_column_number_density_validity > 50; derive(tropospheric_NO2_column_number_density [Pmolec/cm2]); derive(datetime_stop {time}); latitude >= 54.506 [degree_north] ; latitude <= 56.942 [degree_north] ; longitude >= 36.357 [degree_east] ; longitude <= 38.854 [degree_east]; bin_spatial(436, 55.506, 0.001, 495, 37.357, 0.001004); derive(latitude {latitude}); derive(longitude {longitude}); keep(NO2_column_number_density, tropospheric_NO2_column_number_density, stratospheric_NO2_column_number_density, NO2_slant_column_number_density, tropopause_pressure, absorbing_aerosol_index, cloud_fraction, sensor_altitude, sensor_azimuth_angle, sensor_zenith_angle, solar_azimuth_angle, solar_zenith_angle)

Example files: S5P_OFFL_L2__NO2____20190101T082329_20190101T100459_06313_01_010202_20190107T103220.zip
S5P_OFFL_L2__NO2____20190102T080428_20190102T094558_06327_01_010202_20190108T095255.zip

Thank you very much for your time and effort,
Pedro

2

What you are specifying is the number of grid ‘edges’. You always have one more edge than you have grid cells.
The reason you have to specify the edges is so you can specify irregular grids.

HARP does not impose coordinate variables. See also the compatibility section of the HARP documentation. You will have to modify the xarray structure yourself to make certain variables coordinate variables.

For sure. The data that you import is composed of fully native python data structures (with numpy for the data arrays). Have you looked at the examples in the documentation?

3

Thank you so much!

Super helpful!

Yes, I have read that. However, I am following this tutorial by RUS and they apply these transformations and get back an object with lat and lon that can be plotted using matplotlib and cartopy on top of the area of interest, which I find very useful to visually check baby steps. I am doing nothing different but still:

  1. I can’t get the coordinates. Even when I look at the harp product (thanks to your last advice) in python I can’t find the coordinates, can you tell me where to find them? I only see this:

    source product = ‘S5P_OFFL_L2__NO2____20190102T080428_20190102T094558_06327_01_010202_20190108T095255.zip’
    history = “2020-08-05T21:03:01Z [harp-1.11] harp.import_product(’…/data/Moscow/L2__NO2___/S5P_OFFL_L2__NO2____20190102T080428_20190102T094558_06327_01_010202_20190108T095255.zip’,operations=‘tropospheric_NO2_column_number_density_validity > 50; derive(tropospheric_NO2_column_number_density [Pmolec/cm2]); derive(datetime_stop {time}); latitude >= 54.506 [degree_north] ; latitude <= 56.942 [degree_north] ; longitude >= 36.357 [degree_east] ; longitude <= 38.854 [degree_east]; bin_spatial(45, 55.506, 0.01, 51, 37.357, 0.01); derive(latitude {latitude}); derive(longitude {longitude}); keep(NO2_column_number_density, tropospheric_NO2_column_number_density, stratospheric_NO2_column_number_density, NO2_slant_column_number_density, tropopause_pressure, absorbing_aerosol_index, cloud_fraction, sensor_altitude, sensor_azimuth_angle, sensor_zenith_angle, solar_azimuth_angle, solar_zenith_angle)’)”

    double sensor_altitude {time=1, latitude=44, longitude=50} [m]
    double solar_zenith_angle {time=1, latitude=44, longitude=50} [degree]
    double solar_azimuth_angle {time=1, latitude=44, longitude=50} [degree]
    double sensor_zenith_angle {time=1, latitude=44, longitude=50} [degree]
    double sensor_azimuth_angle {time=1, latitude=44, longitude=50} [degree]
    double tropospheric_NO2_column_number_density {time=1, latitude=44, longitude=50} [Pmolec/cm2]
    double NO2_column_number_density {time=1, latitude=44, longitude=50} [mol/m^2]
    double stratospheric_NO2_column_number_density {time=1, latitude=44, longitude=50} [mol/m^2]
    double NO2_slant_column_number_density {time=1, latitude=44, longitude=50} [mol/m^2]
    double cloud_fraction {time=1, latitude=44, longitude=50} []
    double absorbing_aerosol_index {time=1, latitude=44, longitude=50} []
    double tropopause_pressure {time=1, latitude=44, longitude=50} [Pa]

  2. Can you point me to other documentation of plotting harp products with python maybe? I was using xarray because it is really easy with it and helps me to verify if I chose the area of interest correctly.

Thanks for your time :slight_smile:

4

Hello again,

The only way I can see latitude (and longitude) is if I add them into the keep(…) function. However, in this case, I only see 44 values for latitude (50 for longitude), and I believe they are not the center of each pixel… These are the values for latitude I get:

array([55.511, 55.521, 55.531, 55.541, 55.551, 55.561, 55.571, 55.581,
       55.591, 55.601, 55.611, 55.621, 55.631, 55.641, 55.651, 55.661,
       55.671, 55.681, 55.691, 55.701, 55.711, 55.721, 55.731, 55.741,
       55.751, 55.761, 55.771, 55.781, 55.791, 55.801, 55.811, 55.821,
       55.831, 55.841, 55.851, 55.861, 55.871, 55.881, 55.891, 55.901,
       55.911, 55.921, 55.931, 55.941])

Which seem not to be aligned with my desired grid bin_spatial(45, 55.506, 0.01, 51, 37.357, 0.01), but with my former latitude filter (latitude >= 54.506 [degree_north] ; latitude <= 56.942 [degree_north];). So now I am not sure if I got my area of interest and I am unable to plot it, can you help me with this too, please?

This is the full query:

tropospheric_NO2_column_number_density_validity > 50; derive(tropospheric_NO2_column_number_density [Pmolec/cm2]); derive(datetime_stop {time}); latitude >= 54.506 [degree_north] ; latitude <= 56.942 [degree_north] ; longitude >= 36.357 [degree_east] ; longitude <= 38.854 [degree_east]; bin_spatial(45, 55.506, 0.01, 51, 37.357, 0.01); derive(latitude {latitude}); derive(longitude {longitude}); keep(latitude, longitude, NO2_column_number_density, tropospheric_NO2_column_number_density, stratospheric_NO2_column_number_density, NO2_slant_column_number_density, tropopause_pressure, absorbing_aerosol_index, cloud_fraction, sensor_altitude, sensor_azimuth_angle, sensor_zenith_angle, solar_azimuth_angle, solar_zenith_angle)

print(harp_L2_L3)
source product = 'S5P_OFFL_L2__NO2____20190101T082329_20190101T100459_06313_01_010202_20190107T103220.zip'
history = "2020-08-06T12:21:05Z [harp-1.11] harp.import_product('../data/Moscow/L2__NO2___/S5P_OFFL_L2__NO2____20190101T082329_20190101T100459_06313_01_010202_20190107T103220.zip',operations='tropospheric_NO2_column_number_density_validity > 50;     derive(tropospheric_NO2_column_number_density [Pmolec/cm2]);     derive(datetime_stop {time});    latitude >= 54.506 [degree_north] ; latitude <= 56.942 [degree_north] ;    longitude >= 36.357 [degree_east] ; longitude <= 38.854 [degree_east];    bin_spatial(45, 55.506, 0.01, 51, 37.357, 0.01);    derive(latitude {latitude}); derive(longitude {longitude});    keep(latitude, longitude, NO2_column_number_density, tropospheric_NO2_column_number_density, stratospheric_NO2_column_number_density,         NO2_slant_column_number_density, tropopause_pressure, absorbing_aerosol_index, cloud_fraction, sensor_altitude,         sensor_azimuth_angle, sensor_zenith_angle, solar_azimuth_angle, solar_zenith_angle)')"

double sensor_altitude {time=1, latitude=44, longitude=50} [m]
double solar_zenith_angle {time=1, latitude=44, longitude=50} [degree]
double solar_azimuth_angle {time=1, latitude=44, longitude=50} [degree]
double sensor_zenith_angle {time=1, latitude=44, longitude=50} [degree]
double sensor_azimuth_angle {time=1, latitude=44, longitude=50} [degree]
double tropospheric_NO2_column_number_density {time=1, latitude=44, longitude=50} [Pmolec/cm2]
double NO2_column_number_density {time=1, latitude=44, longitude=50} [mol/m^2]
double stratospheric_NO2_column_number_density {time=1, latitude=44, longitude=50} [mol/m^2]
double NO2_slant_column_number_density {time=1, latitude=44, longitude=50} [mol/m^2]
double cloud_fraction {time=1, latitude=44, longitude=50} []
double absorbing_aerosol_index {time=1, latitude=44, longitude=50} []
double tropopause_pressure {time=1, latitude=44, longitude=50} [Pa]
double latitude {latitude=44} [degree_north]
double longitude {longitude=50} [degree_east]

As an example, you can use the following product: S5P_OFFL_L2__NO2____20190101T082329_20190101T100459_06313_01_010202_20190107T103220.zip

Thank you so much,
Pedro H

5

First, you should include the latitude_bounds and longitude_bounds variables in your keep() operation (you want the pixel bounding boxes, not the center points, when you regrid the data).

The size of your final grid is determined by your bin_spatial(45, 55.506, 0.01, 51, 37.357, 0.01) operation. In which you clearly indicated that you wanted a grid with 45 latitude edges (=44 cells) and 51 longitude edges (=50 cells).

6

Dear @ svniemeijer,

Thank you so much! I got what I was looking for. In case you are curious or if you have any suggestions or even if other people are looking for the same, my code can be found in the file ‘mk_raster.py’ in the following repo: https://github.com/pherrusa7/Sentinel5P

Thank you so much!

7

Hi - can you clarify how you were able to get the coordinates in the Xarray dataset? Did you just include them in the keep() variable like the previous user recommended? I’m having the same issue when working through the RUS tutorial.

Thanks!