Cartopy: 10m land feature read incorrectly by basic reader
Description
Depending on the extent of a map, adding the land features will color the entire plot (including the ocean).
Related to #803, #579, #492.
This issue seems to have appeared on the 0.18 version, although I haven't tested it.
The issue is present when plotting with the NearsidePerspective projection, but not PlateCarree (I haven't tried others).
Code to reproduce
import cartopy
import matplotlib.pyplot as plt
nearside = cartopy.crs.NearsidePerspective(central_longitude=0, central_latitude=0)
geo = cartopy.crs.PlateCarree()
ax = plt.axes(projection=nearside)
extent = 7
ax.set_extent([-extent, extent, -extent, extent], crs=geo)
ax.add_feature(cartopy.feature.LAND)
ax.coastlines()
produces
setting the extent to 8 in the above code gives the expected result :
Full environment definition
Operating system
Windows 10
Cartopy version
0.18.0
massoncl
All 23 comments
Which versions of PROJ and GEOS do you have? I am unable to reproduce this.
greglucas
on 12 Jun 2020
I'm using PROJ 7.0.0 and GEOS 3.8.1. I installed everything in a new miniconda environment, with the following commands :
conda create --name test_cartopy python=3.8.3
conda activate test_cartopy
conda install cartopy
Here's the output of conda list, just in case :
brotlipy 0.7.0 py38he774522_1000
bzip2 1.0.8 he774522_0
ca-certificates 2020.4.5.2 hecda079_0 conda-forge
cartopy 0.18.0 py38hde6099d_1 conda-forge
certifi 2020.4.5.2 py38h32f6830_0 conda-forge
cffi 1.14.0 py38h7a1dbc1_0
cfitsio 3.470 he774522_5
chardet 3.0.4 py38_1003
cryptography 2.9.2 py38h7a1dbc1_0
curl 7.69.1 h2a8f88b_0
cycler 0.10.0 py38_0
expat 2.2.9 he025d50_2 conda-forge
freetype 2.10.2 hd328e21_0 conda-forge
freexl 1.0.5 hfa6e2cd_0
gdal 3.1.0 py38h3ba59e7_3 conda-forge
geos 3.8.1 h33f27b4_0
geotiff 1.6.0 hcadc5cd_0 conda-forge
gettext 0.19.8.1 hb01d8f6_1002 conda-forge
glib 2.64.3 he4de6d7_0 conda-forge
hdf4 4.2.13 h712560f_2
hdf5 1.10.6 nompi_ha405e13_100 conda-forge
icc_rt 2019.0.0 h0cc432a_1
icu 64.2 he025d50_1 conda-forge
idna 2.9 py_1
intel-openmp 2020.0 166
jpeg 9d he774522_0 conda-forge
kealib 1.4.13 h3b59ab9_1 conda-forge
kiwisolver 1.2.0 py38h74a9793_0
krb5 1.17.1 hc04afaa_0
libblas 3.8.0 15_mkl conda-forge
libboost 1.67.0 hd9e427e_4
libcblas 3.8.0 15_mkl conda-forge
libcurl 7.69.1 h2a8f88b_0
libffi 3.2.1 h6538335_1007 conda-forge
libgdal 3.1.0 hbc6966c_3 conda-forge
libiconv 1.15 h1df5818_7
libkml 1.3.0 he5f2a48_4
liblapack 3.8.0 15_mkl conda-forge
libnetcdf 4.7.4 nompi_h256d12c_104 conda-forge
libpng 1.6.37 h2a8f88b_0
libpq 12.2 h3235a2c_0
libspatialite 4.3.0a h51df0ed_1038 conda-forge
libssh2 1.9.0 h7a1dbc1_1
libtiff 4.1.0 h56a325e_1
libwebp-base 1.1.0 hfa6e2cd_3 conda-forge
libxml2 2.9.10 h5d81f1c_1 conda-forge
lz4-c 1.9.2 h62dcd97_0
m2w64-expat 2.1.1 2
m2w64-gcc-libgfortran 5.3.0 6
m2w64-gcc-libs 5.3.0 7
m2w64-gcc-libs-core 5.3.0 7
m2w64-gettext 0.19.7 2
m2w64-gmp 6.1.0 2
m2w64-libiconv 1.14 6
m2w64-libwinpthread-git 5.0.0.4634.697f757 2
m2w64-xz 5.2.2 2
matplotlib-base 3.2.1 py38h1626042_0 conda-forge
mkl 2020.0 166
msys2-conda-epoch 20160418 1
numpy 1.18.5 py38h72c728b_0 conda-forge
olefile 0.46 py_0
openjpeg 2.3.1 h57dd2e7_3 conda-forge
openssl 1.1.1g he774522_0
owslib 0.20.0 py_0 conda-forge
pcre 8.44 h6538335_0 conda-forge
pillow 7.1.2 py38hcc1f983_0
pip 20.1.1 py_1 conda-forge
poppler 0.89.0 h0cd1227_0 conda-forge
poppler-data 0.4.9 1 conda-forge
postgresql 12.2 h3235a2c_0
proj 7.0.0 haa36216_3 conda-forge
pycparser 2.20 py_0
pyepsg 0.4.0 py_0
pyopenssl 19.1.0 py38_0
pyparsing 2.4.7 py_0
pyproj 2.6.1.post1 py38h1dd9442_0 conda-forge
pyshp 2.1.0 py_0
pysocks 1.7.1 py38_0
python 3.8.3 he1778fa_0
python-dateutil 2.8.1 py_0
python_abi 3.8 1_cp38 conda-forge
pytz 2020.1 py_0
pyyaml 5.3.1 py38h9de7a3e_0 conda-forge
requests 2.23.0 py38_0
scipy 1.3.2 py38h582fac2_0 conda-forge
setuptools 47.1.1 py38_0
shapely 1.7.0 py38hbf43935_3 conda-forge
six 1.15.0 py_0
sqlite 3.31.1 h2a8f88b_1
tbb 2020.1 he980bc4_0 conda-forge
tiledb 2.0.4 h0b90766_0 conda-forge
tk 8.6.10 hfa6e2cd_0 conda-forge
tornado 6.0.4 py38he774522_1
urllib3 1.25.9 py_0 conda-forge
vc 14.1 h0510ff6_4
vs2015_runtime 14.16.27012 hf0eaf9b_2
wheel 0.34.2 py38_0
win_inet_pton 1.1.0 py38_0
wincertstore 0.2 py38_0
xerces-c 3.2.2 ha925a31_0
xz 5.2.5 h62dcd97_0
yaml 0.2.5 he774522_0 conda-forge
zlib 1.2.11 h62dcd97_4
zstd 1.4.4 ha9fde0e_3
I just ran through some other environment differences between us to check the differences and it looks like the reason this occurs is due to the Shapefile reading differences between pyshp and fiona. I have no idea why they are reading differently, but I would suggest to use fiona if you can it is MUCH faster too.
This environment works for me: conda create --name test_cartopy python=3.8.3 cartopy fiona
greglucas
on 13 Jun 2020
I installed fiona and it solves the issue. Installing fiona with conda downgraded some packages :
The following NEW packages will be DOWNGRADED:
cartopy 0.18.0-py38hde6099d_1 --> 0.18.0-py38h4990f70_0
gdal 3.1.0-py38h3ba59e7_3 --> 3.0.4-py38h3ba59e7_9
geotiff 1.6.0-hcadc5cd_0 --> 1.5.1-h3d29ae3_10
libgdal 3.1.0-hbc6966c_3 --> 3.0.4-h6f60a84_9
poppler 0.89.0-h0cd1227_0 --> 0.87.0-h0cd1227_1
tiledb 2.0.4-h0b90766_0 --> 1.7.7-h0b90766_3
Regardless of whether there's something fishy in the reading, I feel like somewhere we have an issue creating properly closed polygons.
dopplershift
on 14 Jul 2020
We're experiencing the same issue (fatiando/verde#269) and it does only happen with Cartopy 0.18. Don't know if it helps narrow down what might have happened.
leouieda
on 16 Jul 2020
In 0.18, those features started using AutoScaler. With extent=7, that's a width/height of 14 degrees, which will use 10m scaling, and with extent=8, it's 16 degrees, which will use 50m scaling. So you need to force the scale for the best comparison, cartopy.feature.LAND.with_scale('10m').
Then a bisect points to 75c04fbd3d308402f74e13f318648e16eaffba26 or #1210.
QuLogic
on 20 Jul 2020
I think what the difference is is that the original code _always_ created the same type of geometry (i.e., whatever self._geometry_factory = GEOMETRY_FACTORIES.get(shapeType) produced), whereas now it's whatever Shapely decides to create, which I _guess_ is not always the same type.
If we compare the types, the old code always produced MultiPolygon, but the new code produces a Polygon on the last geometry.
QuLogic
on 20 Jul 2020
OK, it's not that, but rather the first MultiPolygon appears to be grouped differently. Haven't quite figured out which one, but it seems to be something outside this region, i.e., the real trigger is using 10m scale, not the small extent or NearSidePerspective.
QuLogic
on 20 Jul 2020
The interior ring of Polygon 4 is missing and is instead the interior ring of Polygon 14. Maybe it's a bug in pyshp's __geo_interface__?
QuLogic
on 20 Jul 2020
I just looked at pyshp's GitHub and the development has stalled it seems, so I'm not sure how likely that would be to get fixed anytime soon.
Unfortunately, if we were to get rid of the BasicReader altogether that would add a hard required dependency on GDAL, which I'm not sure we want to do. However, it seems like a lot of these shapefile releated bugs, my goto solution has been: install Fiona
greglucas
on 20 Jul 2020
So that leaves vendoring pyshp or requiring fiona (and GDAL). I'm not wild about either of those options. I wonder how hard this would be to fix in pyshp.
dopplershift
on 21 Jul 2020
We could override just MultiPolygon conversion. Our old code went through all exterior rings and inserted the interior rings on the ones that actually contained it. The problem is that pyshp places interior rings on the exterior ring that precedes it. However, that is not at all the right polygon in the Natural Earth dataset.
QuLogic
on 21 Jul 2020
Hey all. So the reason for the bug on the part of PyShp, was mostly due to a hack that assumes the ordering of polygon rings, so as to avoid the more expensive hole-in-polygon testing. In my experience, this has rendered fine for most shapefiles, but I agree that the it would be best to follow the shapefile spec more literally to avoid bugs like this. PyShp is still in development, but as the only current maintainer I mostly focus on any major bug fixes, and unfortunately it sometimes takes me a while to get around to it. But following up now with @QuLogic and happy to fix this one.
As a sidenote, I noticed cartopy currently only requires pyshp 1.1.4, but there has since been a major version update to 2.x, with at least some API changes, and a ton of bug fixes. I would recommend updating the requirements to 2.x, especially if you want to be able to avoid the bug that resulted in this issue, once it's taken care of.
karimbahgat
on 5 Aug 2020
So I added a PyShp fix now to correctly assign holes to exteriors based on fast point-in-polygon testing when constructing GeoJSON multipolygons, and tested that it fixed the wrongly assigned hole from this issue (the Black Sea was assigned to the wrong continent). Could someone confirm by using the latest Cartopy and PyShp repos? I'm on Windows and having trouble installing the latest Cartopy version.
karimbahgat
on 14 Aug 2020
Using pip install git+https://github.com/GeospatialPython/pyshp@master, I can confirm this fixes the problem.
QuLogic
on 20 Aug 2020
So then the question becomes, how much effort will it be to update to pyshp 2.x?
dopplershift
on 21 Aug 2020
I'm thinking sometime early next week. I just recently used pyshp on a huge shapefile and caught a small bug in the new code, and added handling of the special case of nested exteriors insides holes etc (eg a small lake on an island inside a lake). Just want to write up some basic tests to confirm with some example shapes and will make a new release this coming week.
karimbahgat
on 22 Aug 2020
Well good to know @karimbahgat, but really I was wondering how much work it would take us to update CartoPy to move to 2.x from 1.x. I just haven't looked to see yet.
dopplershift
on 22 Aug 2020
@dopplershift from my look at cartopy.io.shapereader.py, there really shouldn't be any changes required to switch to pyshp 2.x, the API is still valid and should work.
karimbahgat
on 26 Aug 2020
@QuLogic I discovered some small errors in my original fix, made it much more efficient, and added tests to confirm that it's handling multipolygon holes correctly, incl some edge cases. Would you mind rerunning your tests for this issue using the latest pyshp version on Github? Once you confirm I'll go ahead and push v2.1.1 to PyPI.
karimbahgat
on 26 Aug 2020
I went ahead and pushed the fix for this in pyshp v2.1.1. Also added a PR (#1653) so that Cartopy requires this newest pyshp version, to avoid the rendering bug.
karimbahgat
on 9 Sep 2020
I believe this was fixed upstream, so I'm going to close this now. Feel free to reopen if using the latest pyshp does not solve the issue.
greglucas
on 12 Dec 2020
Related issues
chengzhuzhang
路
8Comments
QuLogic
路
6Comments
maxnoe
路
5Comments
greglucas
路
7Comments
choldgraf
路
7Comments
Most helpful comment
Hey all. So the reason for the bug on the part of PyShp, was mostly due to a hack that assumes the ordering of polygon rings, so as to avoid the more expensive hole-in-polygon testing. In my experience, this has rendered fine for most shapefiles, but I agree that the it would be best to follow the shapefile spec more literally to avoid bugs like this. PyShp is still in development, but as the only current maintainer I mostly focus on any major bug fixes, and unfortunately it sometimes takes me a while to get around to it. But following up now with @QuLogic and happy to fix this one.
As a sidenote, I noticed cartopy currently only requires pyshp 1.1.4, but there has since been a major version update to 2.x, with at least some API changes, and a ton of bug fixes. I would recommend updating the requirements to 2.x, especially if you want to be able to avoid the bug that resulted in this issue, once it's taken care of.