Comments (7)
There is an old PROJ mailing list thread started by me where I discuss this +over
with one of the PROJ developers. My understanding is that +over
should not be used and may be considered deprecated. Instead, it was recommended to me to use +pm=180
or some other equivalent to basically define the projection with the prime meridian not at 0 but instead at 180 degrees. This means your data/polygon values have a contiguous -20 to 20 degrees (pm=180) rather than 160 to -160 (pm=0).
So the longitudes are in the range -180 to 360 technically? Basically going from ~170 to ~220 to keep the pixel offset always positively increasing. Or wait, who is doing the wrapping here? What are the longitudes produced by the AreaDefinition? What projection are you using for mapping/drawing?
from pyresample.
I don't now if this bug is also relevant to the getting full lat/lon arrays. The reason for the bug is that for the boundary calculation in lons/lats the bounding box is calculated in projection coords (geostationary) and then "reprojected" using Proj
(
pyresample/pyresample/geometry.py
Line 2823 in 07382cc
AreaDefinition
.
I also got the feeling that the usage of +over
was not really encouraged but it was the only way I found to get to the "correct" results. In any case its not only a matter of mapping/drawing (this was only to make the problem clear) because if you look at the coordinates and the shapely polygon it is obvious that the boundary is wrong. So if I would do a spatial selection/sjoin I would get wrong results.
While the pm=180 does not mess up the geometry of the polygon the coordinates are shifted and downstream tasks like mapping/selecting are rather difficult and counterintuitive. I think most of the time the recommendation is to split the polygons at the dateline (https://towardsdatascience.com/around-the-world-in-80-lines-crossing-the-antimeridian-with-python-and-shapely-c87c9b6e1513). I don't know what the best way to go here is. I guess this also depends on the usage of the boundary
function.
While a different discussion and maybe special to my usage; but I think it would make sense if boundary
would return a shapely polygon instead of coordinate arrays. That would make it easier to work with (e.g. plotting/ geopandas).
from pyresample.
Some of the work by @ghiggi handles these cases by internally having a Boundary object that contains two shapely polygons for the dateline cases. I don't remember where that implementation lives between future boundary/spherical polygon classes and current classes.
I think I'd prefer if the lon/lat coordinate generation for the bounds did a wrap check manually (especially since bounding coordinates should not have very many elements). So if max(lon) - min(lon)
is greater than 355 and we're not touching the north or south pole, then let's add 360 to the < 0
coordinates. I would at least expect that to work better than what is currently produced and doesn't depend on possibly deprecated PROJ flags. I'm not sure if this would break @ghiggi's checks in the boundary code.
What do you think?
from pyresample.
I am a bit tired right now, but I think the issue here @BENR0 is something else. I guess your Polygon vertices are defined clockwise, while shapely expects to be defined as counterclockwise. Can you try to plot with counterclockwise shapely polygon vertices (with i.e. [, ::-1]) and use the the cartopy Geodetic CRS?
from pyresample.
@djhoese in the future spherical polygon interface I found a way to not need two polygons in the dateline cases :D
I just need to find the time to prepare the PR and implement the test units ...
from pyresample.
@djhoese I agree. The PROJ flags solution, while working or at least giving the results I wanted, seems a little fishy. A manual check would be better I think. From what I grasp I am not sure if it s possible to get away with a single polygon instead of two if it would be crossing the date line but I am very curious what you've come up with @ghiggi.
@ghiggi the counterclockwise solution won't work, at least for the "distorted" shapely polygon display, because the one coordinate that is over the dateline and thus positive is due to the inverse projection using Proj
and shapely does not know anything about the projection/wrapping.
Nevertheless I found out something interesting; the geo plot you see above is made from a geopandas dataframe using hvplot using bokeh which is wrong obviously. I now also did a plot just with cartopy and in that case the polygon looks right. Interestingly this is the case in the clockwise and counterclockwise case even though there is some weird artefact with the background map in the clockwise case. But still cartopy obviously knows how to figure out the dateline crossing.
Spatial join also seem to be correct in the clockwise and counterclockwise cases
So I guess this is a bug with hvplot/bokeh somehow which mislead me.
I don't know if this is still a bug then? I guess it would make sense to change the default to counterclockwise since this is also what shapely expects?
from pyresample.
All pyresample polygon spherical routines expect clockwise-ordered vertices ... and would be an ultra-pain to change.
To what you refer with spatial join? Because the vertices orders defines what is the outside and inside of a polygon.
And cartopy expects counterclockwise vertices. You don't see the weird behavior with clockwise vertices because you just plot the border of the polygon ... but try to fill it and you will see the problem 😄
Did you try to use Geoviews for the plots? You would still rely on bokeh but you could specify the cartopy geodetic crs and solve the issue I guess ....
from pyresample.
Related Issues (20)
- 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
- 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 5
- 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 HOT 1
- RuntimeWarning: invalid value encountered in divide x_2 = (-b__ - np.sqrt(discriminant)) / (2 * a__) HOT 3
- Error in SwathDefinition html representation if lon/lat arrays are dask arrays HOT 5
- Gradient search resampling swath data gives transposed results HOT 4
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.