3

I am drawing custom shape boundary map which focuses on the Pacific Sector of Southern Ocean (160E~180~60W,-60S~-90S). When adding gridlines, the region (180~60W) lacks gridlines.

Example output showing the problem

I have tried many solutions (mostly for rectangular projection), but failed.

import matplotlib.pyplot as plt
import matplotlib.path as mpath
import cartopy
import cartopy.crs as ccrs
import numpy as np
proj = ccrs.Stereographic(central_longitude=228, central_latitude=-70)
fig = plt.figure(figsize=(10,10))
ax = plt.axes([0,0,1,1],projection=proj)
latmin = -90
latmax = -60
lonmin = 150
lonmax = 305
lats = np.linspace(latmax, latmin, latmax - latmin + 1)
lons = np.linspace(lonmin, lonmax, lonmax - lonmin + 1)
vertices = [(lon, latmin) for lon in range(lonmin, lonmax + 1, 1)] + \
    [(lon, latmax) for lon in range(lonmax, lonmin - 1, -1)]
boundary = mpath.Path(vertices)
ax.set_boundary(boundary, transform=ccrs.PlateCarree())
ax.set_extent([170, 320, -90, -60], ccrs.PlateCarree())
ax.add_feature(cartopy.feature.LAND, zorder=1, edgecolor='k')
ax.gridlines(draw_labels=True)
plt.show()
python_user
  • 5,375
  • 2
  • 13
  • 32
Jellywzx
  • 31
  • 3

2 Answers2

2

Gridlines plotting is problematic in this particular case. The plot involves special boundary, the parallels crossing the international dateline -- these may be the causes of the problem. I hope the answer by @Rutger_Kassies is what you need. However, if you still need the original plot extent, here is another approach.

import matplotlib.pyplot as plt
import matplotlib.path as mpath
import cartopy
import cartopy.crs as ccrs
import numpy as np

proj = ccrs.Stereographic(central_longitude=228, central_latitude=-70)
fig = plt.figure(figsize=(10/2, 10/2))
ax = plt.axes([0,0,1,1], projection=proj)

latmin = -90
latmax = -60
lonmin = 150
lonmax = 305
lats = np.linspace(latmax, latmin, latmax - latmin + 1)
lons = np.linspace(lonmin, lonmax, lonmax - lonmin + 1)
vertices = [(lon, latmin) for lon in range(lonmin, lonmax + 1, 1)] + \
    [(lon, latmax) for lon in range(lonmax, lonmin - 1, -1)]

boundary = mpath.Path(vertices)
ax.set_boundary(boundary, transform=ccrs.PlateCarree())

ax.set_extent([170, 320, -90, -60], ccrs.PlateCarree())
ax.add_feature(cartopy.feature.LAND, zorder=1, edgecolor='k')
#ax.add_feature(cartopy.feature.OCEAN, zorder=1, edgecolor='none', alpha=0.3)

# Build gridlines in 2 steps
# First set of gridlines: meridians only
# proper `xlim` is needed to get all the x-labels' visibility set correctly
gl1 = ax.gridlines(draw_labels=True, crs=ccrs.PlateCarree(), \
            xlocs=range(-180, 359, 20), \
            ylocs=[], \
            xlim=[-2847619, 3537282])  # get these from ax.get_xbound() of previous run

# Second set of gridlines: parallels of latitude only. 
# The lines are terminated at the meridian of the dateline!
gl2 = ax.gridlines(draw_labels=True, crs=ccrs.PlateCarree(), \
            y_inline=True, \
            ylocs=range(-80, -50, 10), \
            xlocs=[])

# Draw the missing parts of the parallel lines at 70, 80 deg_S
# Colors are set to `red` and `green` intentionally
lons = range(-180, -50, 5)
ax.plot(lons, -80*np.ones(len(lons)), transform=ccrs.Geodetic(), lw=0.5, color="red", zorder=30)
ax.plot(lons, -70*np.ones(len(lons)), transform=ccrs.Geodetic(), lw=0.5, color="green", zorder=30)

# Print the bounds of the plot
print("Bounds:", ax.get_xbound(), ax.get_ybound())
plt.show()

spole

swatchai
  • 17,400
  • 3
  • 39
  • 58
  • Also check this link: https://github.com/SciTools/cartopy/issues/697#issuecomment-157025400 look at pelson's comment for another approach. – swatchai Apr 28 '23 at 01:38
0

In answer https://stackoverflow.com/a/65690841/4265407 described another method of boundary projection which allows grid to plot correctly.

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

proj = ccrs.Stereographic(central_longitude=228, central_latitude=-70)
fig = plt.figure(figsize=(15,8))
ax = plt.axes([0,0,1,1],projection=proj)

ylim = [-90,-60]
xlim = [150,305]

rect = mpath.Path([[xlim[0], ylim[0]],
                   [xlim[1], ylim[0]],
                   [xlim[1], ylim[1]],
                   [xlim[0], ylim[1]],
                   [xlim[0], ylim[0]],
                   ]).interpolated(40)

proj_to_data   = ccrs.PlateCarree()._as_mpl_transform(ax) - ax.transData
rect_in_target = proj_to_data.transform_path(rect)
ax.set_boundary(rect_in_target)

ax.set_extent([150,305,-90,-58], ccrs.PlateCarree())
ax.add_feature(cartopy.feature.LAND, zorder=1, edgecolor='k')
ax.gridlines(draw_labels=True, x_inline=False, y_inline=False,
             xlocs=range(-180, 180, 10))
plt.show()

map plotting result

Stanislav Ivanov
  • 1,854
  • 1
  • 16
  • 22