Comments (15)
Yes, the order of the vertices is important. A polygon is defined with vertices in a clockwise order around it IIRC.
But the area should never be negative, thanks for reporting this.
from pyresample.
Thanks for the clarification! According to this source, the area of a spherical triangle with angles A,B,C is
R^2*((A+B+C) - PI)
However, the current implementation in get_polygon_area()
is
area += R**2 * e__ - math.pi
which is equivalent to
R^2*(A+B+C) - PI
Maybe that is a problem?
from pyresample.
Ah no, as long as R=1 there is no difference. But for the future it might be good to fix that:
def get_polygon_area(corners):
"""Get the area of the convex area defined by *corners*.
"""
# We assume the earth is spherical !!!
# Should be the radius of the earth at the observed position
R = 1
c1_ = corners[0]
area = 0
for idx in range(1, len(corners) - 1):
b1_ = Arc(c1_, corners[idx])
b2_ = Arc(c1_, corners[idx + 1])
b3_ = Arc(corners[idx], corners[idx + 1])
e__ = (abs(b1_.angle(b2_)) +
abs(b2_.angle(b3_)) +
abs(b3_.angle(b1_)))
area += e__ - math.pi
return R**2 * area
Shall I send a pull request?
from pyresample.
A PR would be great.
Thanks !
from pyresample.
I think I also found the reason for negative polygon areas: In spherical_geometry.Arc.angle()
the angle between two arcs is computed using the dot product of normal vectors ua_
and ub_
:
val = ua_.dot(ub_) / (ua_.norm() * ub_.norm())
if abs(val - 1) < EPSILON:
angle = 0
elif abs(val + 1) < EPSILON:
angle = math.pi
else:
angle = math.acos(val)
If the angle is very small, i.e. the dot product is almost 1, the angle is reset to zero in the subsequent if
block. That causes the sum e__
of angles in the spherical triangle to be slightly less than pi which finally leads to a negative polygon area.
I assume the EPSILON criteria were implemented to prevent dot products like 1.00...1 from triggering math domain errors in math.acos()
? In that case I would propose the following solution which only takes dot products > 1 and < -1 into account:
if 0 <= val - 1 < EPSILON:
angle = 0
elif -EPSILON < val + 1 <= 0:
angle = math.pi
else:
angle = math.acos(val)
from pyresample.
Well actually everything outside of [-1, 1] would trigger an error, so maybe we could simply write
if val > 1: angle = 0
etc... ?
from pyresample.
That's true. But it would mask a possible bug in the dot product computation: Even if the dot product yielded a completely unrealistic value, the angles would still be valid.
from pyresample.
Good point.
Also, I remember now the use of EPSILON: it's also to be able to detect colinearity, so smaller angles are snapped to 0...
from pyresample.
I see. What about adding a snap
argument to the angle()
method,
if snap:
if abs(val - 1) < EPSILON:
angle = 0
elif abs(val + 1) < EPSILON:
angle = math.pi
else:
angle = math.acos(val)
else:
if 0 <= val - 1 < EPSILON:
angle = 0
elif -EPSILON < val + 1 <= 0:
angle = math.pi
else:
angle = math.acos(val)
and calling it with snap=False
in get_polygon_area()
?
from pyresample.
Sounds good 👍
from pyresample.
The last two tests on overlap rate would then fail:
======================================================================
FAIL: Test how much two areas overlap.
----------------------------------------------------------------------
Traceback (most recent call last):
File "/home/sfinkens/code/pyresample/pyresample/test/test_spherical_geometry.py", line 154, in test_overlap_rate
self.assertAlmostEqual(area1.overlap_rate(area2), 0.5, 2)
AssertionError: 0.5086952608806948 != 0.5 within 2 places
======================================================================
FAIL: Test how much two areas overlap.
----------------------------------------------------------------------
Traceback (most recent call last):
File "/home/sfinkens/code/pyresample/pyresample/test/test_spherical_geometry.py", line 155, in test_overlap_rate
self.assertAlmostEqual(area2.overlap_rate(area1), 0.068, 3)
AssertionError: 0.06850857188651323 != 0.068 within 3 places
Although the computed values are very close. Original values are 0.499837070923 and 0.0679026012212. Rounding to (1,2) instead of (2,3) digits would resolve the issue, but I don't know whether this is correct.
from pyresample.
Looking at the test code, I don't think there is a reason 0.499837070923 should be more correct than 0.5086952608806948. And since you removed the snapping, I think the latter should be preferred.
What about changing the test values to 0.509 and 0.0685 ?
from pyresample.
Ok!
from pyresample.
Merged, thanks for bringing this up !
from pyresample.
Happy to help. Thanks for the great work you have done already!
from pyresample.
Related Issues (20)
- `get_neighbour_info` slows down significantly when working with large target rasters using many segments HOT 3
- `gradient_search` fails when resampling Himawari data HOT 7
- EWA resampling in 1.27 slows down four times than 1.26.1 HOT 69
- Refactor boundary creation logic
- area definition for a rotated pole coordinate system HOT 24
- Index Error when calling `boundary` with non full disk geos ara HOT 5
- Bug where boundary coordinates wrap around 180 degrees when geostationary AreaDefinition extends beyond 180 degrees HOT 7
- Upgrade to Cython 3.0 and check annotations
- Compatibility with libproj v9.3 HOT 23
- How should this warning be addressed? HOT 2
- Catch or fix pyproj UserWarning when loading an AreaDefinition from a netCDF/CF file HOT 1
- Make a Swath definition subclass for interpolated data HOT 4
- The finding of the optimal radius of influence makes assumption on the ordering of the dimensions in the longitude array HOT 4
- Wrong coordinates returned by `AreaDefintion.get_lonlats` for some projections in out-of-Earth locations HOT 6
- dump DynamicAreaDefinition to YAML
- errors in area definition should not be silently ignored HOT 1
- Resampling GOES mesoscale data to my area gives blank data HOT 10
- Remove Configobj as Dependency Due To Security Vulnerability and No Longer Maintaned HOT 3
- YAML area configuration does not allow to specify the dtype and dtype is ignored in equality comparisons HOT 6
- "resample_nearest" GEO swath to eqc has tiny periodic "glitch" in results
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from pyresample.