# Problems with Harp on Pyhton
**URL:** https://forum.atmospherictoolbox.org/t/problems-with-harp-on-pyhton/195
**Category:** HARP
**Tags:** sentinel-5p, python
**Created:** [August 4, 2020, 5:46pm UTC](https://forum.atmospherictoolbox.org/t/problems-with-harp-on-pyhton/195 "2020-08-04T17:46:56Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@pherrusa7](https://forum.atmospherictoolbox.org/u/pherrusa7)
#### Post date: [August 4, 2020, 5:46pm UTC](https://forum.atmospherictoolbox.org/t/problems-with-harp-on-pyhton/195/1 "2020-08-04T17:46:56Z")
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
---
### Author:  [@sander.niemeijer](https://forum.atmospherictoolbox.org/u/sander.niemeijer)
#### Post date: [August 5, 2020, 11:02am UTC](https://forum.atmospherictoolbox.org/t/problems-with-harp-on-pyhton/195/2 "2020-08-05T11:02:00Z")
> [@pherrusa7](#):
>
> - 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.
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.
> [@pherrusa7](#):
>
> - 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?
HARP does not impose coordinate variables. See also the [compatibility section of the HARP documentation](http://stcorp.github.io/harp/doc/html/conventions/compatibility.html). You will have to modify the xarray structure yourself to make certain variables coordinate variables.
> [@pherrusa7](#):
>
> - 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?
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](http://stcorp.github.io/harp/doc/html/python.html#examples)?
---
### Author:  [@pherrusa7](https://forum.atmospherictoolbox.org/u/pherrusa7)
#### Post date: [August 5, 2020, 10:11pm UTC](https://forum.atmospherictoolbox.org/t/problems-with-harp-on-pyhton/195/3 "2020-08-05T22:11:09Z")
Thank you so much!
> [@svniemeijer](#):
>
> 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.
Super helpful!
> [@svniemeijer](#):
>
> HARP does not impose coordinate variables. See also the [compatibility section of the HARP documentation](http://stcorp.github.io/harp/doc/html/conventions/compatibility.html). You will have to modify the xarray structure yourself to make certain variables coordinate variables.
Yes, I have read that. However, I am following this [tutorial](https://www.youtube.com/watch?v=CE6BeLPORIE) 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:
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 🙂
---
### Author:  [@pherrusa7](https://forum.atmospherictoolbox.org/u/pherrusa7)
#### Post date: [August 6, 2020, 12:40pm UTC](https://forum.atmospherictoolbox.org/t/problems-with-harp-on-pyhton/195/4 "2020-08-06T12:40:41Z")
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
---
### Author:  [@sander.niemeijer](https://forum.atmospherictoolbox.org/u/sander.niemeijer)
#### Post date: [August 7, 2020, 9:26am UTC](https://forum.atmospherictoolbox.org/t/problems-with-harp-on-pyhton/195/5 "2020-08-07T09:26:26Z")
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).
---
### Author:  [@pherrusa7](https://forum.atmospherictoolbox.org/u/pherrusa7)
#### Post date: [August 17, 2020, 11:38am UTC](https://forum.atmospherictoolbox.org/t/problems-with-harp-on-pyhton/195/6 "2020-08-17T11:38:56Z")
Dear @ [svniemeijer](https://forum.atmospherictoolbox.org/u/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](https://github.com/pherrusa7/Sentinel5P)
Thank you so much!
---
### Author:  [@priyapatel1](https://forum.atmospherictoolbox.org/u/priyapatel1)
#### Post date: [May 6, 2021, 6:25pm UTC](https://forum.atmospherictoolbox.org/t/problems-with-harp-on-pyhton/195/7 "2021-05-06T18:25:15Z")
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!