Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

No warning if extents past limits are clipped #1730

Open
kdpenner opened this issue Feb 3, 2021 · 6 comments
Open

No warning if extents past limits are clipped #1730

kdpenner opened this issue Feb 3, 2021 · 6 comments

Comments

@kdpenner
Copy link
Contributor

kdpenner commented Feb 3, 2021

Description

TransverseMercator has two properties bounding map units:

cartopy/lib/cartopy/crs.py

Lines 856 to 862 in 2148e13

@property
def x_limits(self):
return (-2e7, 2e7)
@property
def y_limits(self):
return (-1e7, 1e7)

I'd like to see past these boundaries, but set_xlim and set_ylim don't work as expected.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

tmerc_crs_plt = ccrs.TransverseMercator(approx=False)
region_proj = [-5e6, 5e6, -5e6, 2e7]
x1, x2, y1, y2 = region_proj

fig = plt.figure(figsize=(20, 7))
ax3 = fig.add_subplot(1, 2, 1, projection=tmerc_crs_plt)
ax3.set_extent(region_proj, crs=tmerc_crs_plt)
ax3.coastlines()

ax5 = fig.add_subplot(1, 2, 2, projection=tmerc_crs_plt)
ax5.set_xlim(x1, x2)
ax5.set_ylim(y1, y2)
ax5.coastlines()

plt.show()

tmerc

On the left, the top is clipped to 1e7 without warning. On the right, setting 2e7 creates whitespace instead of moving the spine and plotting the additional data.

It seems like these properties are used all over---what are the consequences of changing them?

@blazing216
Copy link
Contributor

The x_limits and y_limits are used to define the boundary. Everything outside the boundary is clipped.

cartopy/lib/cartopy/crs.py

Lines 848 to 854 in 7c75f5a

@property
def boundary(self):
x0, x1 = self.x_limits
y0, y1 = self.y_limits
return sgeom.LinearRing([(x0, y0), (x0, y1),
(x1, y1), (x1, y0),
(x0, y0)])

I think you have to change the value hard-coded in the x_limits and y_limits if you want to go beyond the limits. The largest I can go is x_limits to (-1.671e7, 1.671e7) and y_limits to (-2e7, 2e7). Beyond that, the map is blank. It gives you something similar to https://en.wikipedia.org/wiki/File:MercTranSph.png
clon_0 0

@dopplershift
Copy link
Contributor

There's definitely a sensible max for these coords, based on the PROJ docs. We could probably do better and make them something like y = +/- scale_factor * radius * np.pi/2 and x = +/- scale_factor * radius * const. I tried to come up with a value for const, but it seems like it depends on how small and increment of points I use, so it doesn't seem to be converging.

@QuLogic
Copy link
Member

QuLogic commented Jul 19, 2021

The defined area from the docs is Global, with full accuracy within 3900 km of the central meridian, though it's only 0.1mm at 7000km. In fact, before Transverse Mercator was replaced with etmerc (or, approx=True right now), there was much more severe distortion as you went out further.

Vertically, Mercator extends to infinity at the pole, but Transverse Mercator is the opposite, and we need to pick a cut off point for x.

@kdpenner kdpenner changed the title Unable to set x and y limits outside CRS's properties No warning if extents past limits are clipped Jul 20, 2021
@kdpenner
Copy link
Contributor Author

I changed the title of the PR to better reflect the issue. I expected to see the extent I passed or a warning, not a silent clip to the hard coded limits. In part because I didn't know they existed.

@dopplershift I changed the limits as part of #1739

@QuLogic where do the docs mention the chosen area? the x and y limits as they are now---x = (-2e7, 2e7), y = (-1e7, 1e7)---include regions that don't have an inverse transformation.

@blazing216
Copy link
Contributor

blazing216 commented Jul 21, 2021

Instead of setting a cut off for x and y, maybe it is better to set the cut off for longitude and latitude.

The reason I proposed that is because some projections, e.g. Transverse Mercator, do not work properly beyond certain limits. They may project some points that should be outside the x, y limits inside. Hence, those points cannot be cut off properly. They are projected incorrectly because you transform those projected points back, the resulting longitude and latitude are not the same as the original ones.

I saw there are some latitude limits set for Mercator. I guess it is used to cut off the geometries before projection?

You may see an example when plotting coastlines with TransverseMercator(central_longitude=-40.9). The line across South America is from a wrong projection of a segment from Africa.

If I understand it correctly, _project_segment uses a clear way to split a line segment at the 'timeline'. However, looks like it assumes the projections work properly for all the longitude and latitudes, which is not true. Maybe adding a check of the validity of the projected points can be helpful.

else:
t_current = t_max
p_current = p_max
state = get_state(p_current, gp_domain, handle)
if state == POINT_IN:
lines.new_line()

clon_-40 9

@kdpenner
Copy link
Contributor Author

@blazing216 Transverse Mercator has holes which are centered on the equator and (-90-central_longitude, 90+central_longitude). A simple cutoff will not work. #1739 fixes your example, but my guess is that a solution in trace.pyx will be better than mine.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants