joernu76
commented
5 years ago Cross-link to related issue of MSS project at https://bitbucket.org/wxmetvis/mss/issues/420/cartopy-bounding-box-setting-with
Open joernu76 opened 5 years ago
joernu76
commented
5 years ago Cross-link to related issue of MSS project at https://bitbucket.org/wxmetvis/mss/issues/420/cartopy-bounding-box-setting-with
kdpenner
commented
3 years ago In digging around to solve another set_extent issue I think I've found the culprit for this one. There seems to be a numerical problem with cartopy.trace.project_linear. Comparing cartopy's performance with pyproj's:
import cartopy
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from pyproj import CRS, Transformer
from pyproj.crs.coordinate_operation import StereographicConversion
from pyproj.crs import ProjectedCRS
from pyproj import Geod
from shapely.geometry import LineString
import numpy as np
lonlat_crs = CRS('EPSG:4326')
stereo = StereographicConversion(latitude_natural_origin=90)
stereo_crs = ProjectedCRS(conversion=stereo)
stereo_transform_func = Transformer.from_crs(
lonlat_crs, stereo_crs, always_xy=True).transform
lonlat_transform_func = Transformer.from_crs(
stereo_crs, lonlat_crs, always_xy=True).transform
region = [45, -135, 0, 0]
x1, x2, y1, y2 = region
bbox = [[x1, y1], [x2, y1], [x2, y2], [x1, y2], [x1, y1]]
g = Geod(ellps="WGS84")
pts_pyproj = []
for i, pts in enumerate(bbox[:-1]):
pts_lonlat = np.array(g.npts(*pts, *bbox[i+1], 7))
pts_lonlat = np.insert(pts_lonlat, 0, pts, axis=0)
xs, ys = stereo_transform_func(*pts_lonlat.T)
pts_pyproj.extend(list(zip(xs, ys)))
pts_pyproj = np.array(pts_pyproj).T
lons_pyproj, lats_pyproj = lonlat_transform_func(*pts_pyproj)
lonlat_crs_plt = ccrs.PlateCarree()
stereo_crs_plt = ccrs.Stereographic(central_longitude=0, central_latitude=90)
xs, ys = cartopy.trace.project_linear(LineString(bbox), lonlat_crs_plt,
stereo_crs_plt)[0].xy
lonlatz_cartopy = lonlat_crs_plt.transform_points(stereo_crs_plt, np.array(xs),
np.array(ys))
lons_cartopy, lats_cartopy = lonlatz_cartopy[:, 0:2].T
fig = plt.figure(figsize=(10, 4))
ax1 = fig.add_subplot(1, 2, 1)
ax1.scatter(*pts_pyproj, color="tab:blue", label="pyproj")
ax1.scatter(xs, ys, color="tab:orange", label="cartopy")
ax1.set_title("extent -> projected")
ax2 = fig.add_subplot(1, 2, 2)
ax2.scatter(lons_pyproj, lats_pyproj, color="tab:blue", label="pyproj")
ax2.scatter(lons_cartopy, lats_cartopy, color="tab:orange", label="cartopy")
ax2.set_title("extent -> projected -> back to lon lat")
ax1.legend()
ax2.legend()
plt.savefig("stereo.png", bbox_inches="tight")
plt.close()produces this:
the extent in projected space is set by the bounds of the geometries in the plot on the left.
kdpenner
commented
3 years ago I take back my previous comment. set_extent behaves differently between a CRS and a geodetic. Instead of passing ccrs.PlateCarree() in set_extent, you need to pass ccrs.PlateCarree().as_geodetic() or omit it.
joernu76
commented
3 years ago Ha! I can confirm that this works. The behaviour and probably the code has changed in the last 18 month as also my unmodified testcase looks "better" than before.
To make the two figures look identical, the as_geodetic() has to be added. This works as expected:
import os
import matplotlib.pyplot as plt
from netCDF4 import Dataset as netcdf_dataset
import numpy as np
from cartopy import config
import cartopy.crs as ccrs
# get the path of the file. It can be found in the repo data directory.
fname = os.path.join(config["repo_data_dir"],
'netcdf', 'HadISST1_SST_update.nc')
dataset = netcdf_dataset(fname)
sst = dataset.variables['sst'][0, :, :]
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]
prj = ccrs.Stereographic(central_longitude=0, central_latitude=90)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=prj)
print("default", ax.get_extent())
ax.set_extent([45, -135, 0, 0], ccrs.PlateCarree().as_geodetic()) # <-- changed here
print("set1", ax.get_extent())
ax.contour(lons, lats, sst, 60, transform=ccrs.PlateCarree())
ax.coastlines()
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=prj)
print("default", ax.get_extent())
coords = prj.transform_points(
ccrs.PlateCarree(), np.asarray([45, -135]), np.asarray([0, 0]))
ax.set_extent([coords[0, 0], coords[1, 0], coords[0, 1], coords[1, 1]], prj)
print("set2", ax.get_extent())
ax.contour(lons, lats, sst, 60, transform=ccrs.PlateCarree())
ax.coastlines()
plt.show()Issue can be closed IMO.
kdpenner
commented
3 years ago The solution here has to do with the choices discussed in #131 and #645, so it might be useful to keep open.
Description
Using set_extent to determine the bounding box of the plot gives strange results for Stereographic projections ad extents given in longitude/latitude coordinates.
The example below demonstrates the behaviour. The extent of the projection shall go from -45E/0N in the lower left cornert to 135E/0N in the upper right corner. The first figure has a correct upper right corner, but in my case the lower left crner is located at -45E/-20N. For the second figure, I compute the desired coordinates of the bound using the projection and apply them in projectio-space.
Expected behaviour: Identical plots in both figures.
Code to reproduce
Traceback
not applicable
Operating system
Linux
Cartopy version
0.16.0