IllegalArgumentException Using Cartopy to Map Upper-air model data

[spurpuraUNCC]

Description

I am using Cartopy to create upper-air atmospheric maps. I have been following issue #1551 and #879, but I am still have the issue. The code runs and I get the map just fine, but it takes like 20 minutes to run the code. I get about 1000 of these error messages when it is running: ERROR:shapely.geos:IllegalArgumentException: Argument must be Polygonal or LinearRing

Code to reproduce

#Modules:
from datetime import datetime
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.units import units
import numpy as np
from scipy.ndimage import gaussian_filter
import xarray as xr
from netCDF4 import Dataset, num2date
from siphon.ncss import NCSS
import pyproj

#Mapping 250 hPa height, wind speed

##Enter year,month,day,hour
dt=datetime(2019, 2, 20, 18)

print(dt)

##will pull out the dataset and use subset to set up requesting a subset of the data
##'We can then use the ncss object to create a new query object, which facilitates asking for data from the server.'
##'We can look at the ncss.variables object to see what variables are available from the dataset:'
ncss = NCSS(f'https://www.ncei.noaa.gov/thredds/ncss/model-rap130anl-old/{dt:%Y%m}/'
f'{dt:%Y%m%d}/rap_130_{dt:%Y%m%d_%H}00_000.grb2')
query = ncss.query()

print('requesting data from thredds server...')
print()

##make query to get data from the lat and long

query.lonlat_box(north=60, south=18, east=300, west=225)
query.all_times()
query.add_lonlat()
query.accept('netcdf').add_lonlat()
query.variables('Geopotential_height_isobaric','u-component_of_wind_isobaric',
'v-component_of_wind_isobaric')

##'We now request data from the server using this query. The NCSS class handles parsing this NetCDF data

data = ncss.get_data(query)

Extract data and assign units

hght = 0
uwnd = 0
vwnd = 0

#extract coordinate data for plotting
lat = data.variables['lat'][:]
lon = data.variables['lon'][:]
lev = 0

#extract data and assign units

hght = gaussian_filter(data.variables['Geopotential_height_isobaric'][0],

sigma=1.0) * units.meter
uwnd = gaussian_filter(data.variables['u-component_of_wind_isobaric'][0], sigma=1.0) * units('m/s')
vwnd = gaussian_filter(data.variables['v-component_of_wind_isobaric'][0], sigma=1.0) * units('m/s')

#extract coordinate data for plotting
lat = data.variables['lat'][:]

lon = data.variables['lon'][:]
lev = data.variables['isobaric'][:]

time = data.variables['time']
vtime = num2date(time[0], units=time.units)

#Calculate dx and dy for calculations 
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)

#specify 500 hPa data
ilev250 = 0
hght_250 = 0
uwnd_250 = 0
vwnd_250 = 0

#specifcy 250 hPa data

ilev250 = np.where(lev == 25000)[0][0]
hght_250 = hght[ilev250]
uwnd_250 = uwnd[ilev250]
vwnd_250 = vwnd[ilev250]

#use MetPy to calculate the wind speed for colorful plot, change units to knots from m/s

sped_250 = mpcalc.wind_speed(uwnd_250, vwnd_250).to('kt')

#plot projection. The look you want for the view, LambertConformal for mid-latitude view

plotproj = ccrs.LambertConformal(central_longitude=-100., central_latitude=40.,
standard_parallels=[30, 60])

datacrs = ccrs.PlateCarree()

def create_map_background():
fig=plt.figure(figsize=(14, 12))
ax=plt.subplot(111, projection=plotproj)
ax.set_extent([-125, -73, 25, 50],ccrs.PlateCarree())
ax.coastlines('50m', linewidth=0.75)
ax.add_feature(cfeature.STATES, linewidth=0.5)
return fig, ax

fig, ax = create_map_background()

#plot 250 hPa colorful wind speed in knots
clevs_250_sped = np.arange(50, 250, 20)
cf = ax.contourf(lon, lat, sped_250, clevs_250_sped, cmap=plt.cm.BuPu,
transform=datacrs)
plt.colorbar(cf, orientation='horizontal', pad=0, aspect=50)

clevs_250 = np.arange(0, 11000, 120)
cs = ax.contour(lon, lat, hght_250, clevs_250, colors='black',
transform=datacrs)
plt.clabel(cs, fmt='%d’)

#make some nice titles for the plot (one right, one left)
plt.title('250-hPa RAP Geopotential Heights (m) and Wind Speed (kt)', loc='left')
plt.title('Valid Time: {}'.format(vtime), loc='right')

#adjust image and show

plt.subplots_adjust(bottom=0, top=1)
plt.savefig('250_'+str(dt)+'_hght_wind.png',dpi=150)

Traceback

ERROR:shapely.geos:IllegalArgumentException: Argument must be Polygonal or LinearRing

Full environment definition

### Operating system ### Cartopy version ### conda list ``` ``` ### pip list ``` ```

I have reinstalled:

pip install Cython Cython==0.29.21 pip install --no-binary :all: shapely Shapley==1.7.1 pip install pyshp pyshp==2.1.2 pip install six six==1.15.0

brew install proj brew install geos

pip install Cartopy Cartopy==2.1.2

[ClaudeDRT]

Thankyou for the input @spurpuraUNCC ! For anyone overlooking this I had the same issue even when completely reinstalling anaconda and first uninstalled shapely and cartopy and then used: pip install shapely conda install -c conda-forge cartopy

the mimimum way of getting the same ERROR is by running the following (no data required):

import cartopy
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_axes([0,0,1,1], projection=ccrs.LambertAzimuthalEqualArea())
ax.set_extent([10,20,10,20], crs=ccrs.PlateCarree())

(IllegalArgumentException: Argument must be Polygonal or LinearRing)

HOWEVER this wont produce the error and cartopy was totally happy with it:

import cartopy
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_axes([0,0,1,1], projection=ccrs.PlateCarree())
ax.set_extent([10,20,10,20], crs=ccrs.PlateCarree())

The same goes for all map projections that are nice rectangular shape and have rectangular griddings. I realised that this is because you can define the edges of the map in these cases just by using [min_lon, max_lon, min_lat, max_lat] because that's enough information to define the four corners.

I think therefore in a general sense the error message means that there isn't enough information in the coordinates that you're providing for cartopy to know what to plot. This is probably because proj/shapely isn't properly installed and linked (a common problem when using brew to install proj/shapely for the latest Mac OSX versions) and so shapely doesn't really know what to do when making polygon objects and so cartopy doesn't know how to make polygons.

What I would definitely recommend is, in line with the dependencies of cartopy 0.18 (the latest version), to make cartopy work (using python 3.8 because I dont think cartopy works with 3.9 yet because of this https://scitools.org.uk/cartopy/docs/latest/whats_new.html#) do the following:

conda create -n cartopy-env python=3.8.5
source activate cartopy-env
conda uninstall -c conda-forge cartopy
conda uninstall -c conda-forge proj shapely
conda install -c conda-forge cartopy==0.18
conda install jupyterlab
cd /Users/yourusername/to/where/your/.ipynb/files/are
jupyter lab

This will install cartopy, proj, geos, shapely and other dependencies (A full list of what is installed is here: c-ares, cartopy, cycler, freetype, geos, jpeg, kiwisolver, krb5, lcms2, libblas, libcblas, libcurl, libedit, libev, libgfortran, libgfortran5, liblapack, libnghttp2, libopenblas, libpng, libssh2, libtiff, libwebp-base, llvm-openmp, lz4-c, matplotlib-base, numpy, olefile, pillow, proj, pyparsing, pyshp, python-dateutil, scipy, shapely, six, tornado and zstd).

The REQUIRED DEPENDENCIES of cartopy are: GEOS, NumPy, Cython, Shapely, proj, pyshp and six, so UNINSTALL THESE and RUN conda install -c conda-forge cartopy==0.18 if you have problems.

Also, just in case, make sure you DON'T have any lines like: export GEOS_DIR=/usr/local/Cellar/geos/3.8.1/ for any of the required dependencies in you .bashrc, .bash_profile or .zshrc file (which you can access using nano ~/.bashrc or nano ~/.zshrc) and if you do then probably delete them (unless you know exactly what you're doing with them) and open a new terminal window (and activate your cartopy environment) and continue before trying to use cartopy.

I made sure that this method works with both cartopy 0.17 and cartopy 0.18, as well as using methods such as ax.set_extent(), ax.gridlines(), ax.coastlines() and ax.add_feature().

I hope this helps someone:)))

[dopplershift]

I really hate this bug.

@ClaudeDRT Thanks for trying to post code to reproduce. Unfortunately, it doesn't crash for me, and honestly I wouldn't expect it to unless you at least plotted some kind of map or something (otherwise, there's nothing that's going to trigger the offending call into Shapely)

@spurpuraUNCC Sorry you've run into this. The only (completely unsatisfying) work-arounds I've seen are to adjust the map extent or tweak contour intervals/levels so that you don't produce the invalid shapes.

[ClaudeDRT]

I really hate this bug.

@ClaudeDRT Thanks for trying to post code to reproduce. Unfortunately, it doesn't crash for me, and honestly I wouldn't expect it to unless you at least plotted some kind of map or something (otherwise, there's nothing that's going to trigger the offending call into Shapely)

@spurpuraUNCC Sorry you've run into this. The only (completely unsatisfying) work-arounds I've seen are to adjust the map extent or tweak contour intervals/levels so that you don't produce the invalid shapes.

@dopplershift I think the 'ERROR' will still actually run without crashing, but will show the 'ERROR' message (which is more like a warning message but not quite a warning because it can't be ignored like a warning can:£) 8 times, unless perhaps anaconda has updated or it only shows in jupyter lab?

[dopplershift]

@ClaudeDRT I can't even make that code fail in jupyter lab with Python 3.9 and the whole stack from conda-forge.

[ClaudeDRT]

@ClaudeDRT I can't even make that code fail in jupyter lab with Python 3.9 and the whole stack from conda-forge.

Apologies I think I must've done pip install shapely before conda install -c conda-forge cartopy for the exception to show.

I just created a new python 3.9 environment, then ran pip install shapely netCDF4, then ran conda install -c conda-forge cartopy and pip install jupyterlab and then ran the exception code and it came up with the following:

However, if I create a new python 3.9 environment and just run conda install -c conda-forge cartopy and then run pip install jupyterlab I can rerun the same code without an exception. I would I'd be cautious in general using python 3.9 instead of python 3.8 with cartopy 0.18 because I don't think it's supported yet but I realise that you're a contributer to cartopy so probably know more than I do @dopplershift ^^).

[dopplershift]

@ClaudeDRT I'm not sure why you don't think CartoPy 0.18 is supported on Python 3.9--the packages are there for conda-forge and work fine.

Honestly, if you're pip installing Shapely but installing CartoPy from conda-forge, I'm utterly amazed you're not seeing a crash. Really, you should just be doing :conda install -c conda-forge cartopy jupyterlab. I'm not sure why you'd need pip for jupyterlab.