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 18 forks source link

Wrong result of area_covers_point #272

Closed zxdawn closed 2 years ago

zxdawn commented 2 years ago

Sometimes the swath covers the specific point, but it raises no point inside error. Here's an example:

Example

import harp
import matplotlib.pyplot as plt

filename = '../data/tropomi/S5P_NRTI_L2__NO2____20220612T215110_20220612T215610_24168_02_020301_20220612T234455.nc'

operations = ";".join([
    "latitude>60",
    # "area_covers_point(100, 55)",
])

product = harp.import_product(filename, operations)

fig, ax = plt.subplots()

ax.scatter(product.longitude.data, product.latitude.data, c=product.NO2_column_number_density.data)
ax.scatter(100, 70, c='r')

image

area_covers_point

operations = ";".join([
    "latitude>60",
    "area_covers_point(100, 70)",
])

product = harp.import_product(filename, operations)

Error:

NoDataError                               Traceback (most recent call last)
/home/xin/Documents/github/S5P-LNO2-Notebook/analysis/lightning_in_swaths.ipynb Cell 9' in <module>
      [3](vscode-notebook-cell:/home/xin/Documents/github/S5P-LNO2-Notebook/analysis/lightning_in_swaths.ipynb#ch0000030?line=2) filename = '../data/tropomi/S5P_NRTI_L2__NO2____20220612T215110_20220612T215610_24168_02_020301_20220612T234455.nc'
      [6](vscode-notebook-cell:/home/xin/Documents/github/S5P-LNO2-Notebook/analysis/lightning_in_swaths.ipynb#ch0000030?line=5) operations = ";".join([
      [7](vscode-notebook-cell:/home/xin/Documents/github/S5P-LNO2-Notebook/analysis/lightning_in_swaths.ipynb#ch0000030?line=6)     "latitude>60",
      [8](vscode-notebook-cell:/home/xin/Documents/github/S5P-LNO2-Notebook/analysis/lightning_in_swaths.ipynb#ch0000030?line=7)     "area_covers_point(100, 70)",
      [9](vscode-notebook-cell:/home/xin/Documents/github/S5P-LNO2-Notebook/analysis/lightning_in_swaths.ipynb#ch0000030?line=8) ])
---> [11](vscode-notebook-cell:/home/xin/Documents/github/S5P-LNO2-Notebook/analysis/lightning_in_swaths.ipynb#ch0000030?line=10) product = harp.import_product(filename, operations)

File ~/miniconda3/lib/python3.9/site-packages/harp/_harppy.py:1305, in import_product(filename, operations, options, reduce_operations, post_operations)
   [1302](file:///home/xin/miniconda3/lib/python3.9/site-packages/harp/_harppy.py?line=1301) try:
   [1303](file:///home/xin/miniconda3/lib/python3.9/site-packages/harp/_harppy.py?line=1302)     # Raise an exception if the imported C product contains no variables, or variables without data.
   [1304](file:///home/xin/miniconda3/lib/python3.9/site-packages/harp/_harppy.py?line=1303)     if _lib.harp_product_is_empty(c_product_ptr[0]) == 1:
-> [1305](file:///home/xin/miniconda3/lib/python3.9/site-packages/harp/_harppy.py?line=1304)         raise NoDataError()
   [1307](file:///home/xin/miniconda3/lib/python3.9/site-packages/harp/_harppy.py?line=1306)     # Convert the C product into its Python representation.
   [1308](file:///home/xin/miniconda3/lib/python3.9/site-packages/harp/_harppy.py?line=1307)     product = _import_product(c_product_ptr[0])

NoDataError: product contains no variables, or variables without data

Own function

I asked a similar question on stackoverflow and the solution works well:

from scipy.optimize import linprog

def in_hull(points, x):
     # https://stackoverflow.com/a/43564754/2912349
     n_points = len(points)
     n_dim = len(x)
     c = np.zeros(n_points)
     A = np.r_[points.T,np.ones((1,n_points))]
     b = np.r_[x, np.ones(1)]
     lp = linprog(c, A_eq=A, b_eq=b)
     return lp.success

points = np.c_[product.longitude.data, product.latitude.data]
in_hull(points, [100, 70])  # returns True
svniemeijer commented 2 years ago

You are mixing latitude and longitude. Please check the documentation. You need to use area_covers_point(latitude [unit], longitude [unit]).