Unable to plot lakes or provinces

[alexandrefierro]

Greetings:

I am able to plot rivers, boundaries, oceans, landuse etc but for some unknown reasons lakes and provinces always return this gibberish below (with res=10, 50, 110 m). Any idea how to bypass this would be very much appreciated. I have tested nearly all versions of cartopy >= 0.17 but in vain.

Traceback (most recent call last): File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/matplotlib/backends/backend_qt5.py", line 506, in _draw_idle self.draw() File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/matplotlib/backends/backend_agg.py", line 388, in draw self.figure.draw(self.renderer) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/matplotlib/artist.py", line 38, in draw_wrapper return draw(artist, renderer, *args, kwargs) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/matplotlib/figure.py", line 1709, in draw renderer, self, artists, self.suppressComposite) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/matplotlib/image.py", line 135, in _draw_list_compositing_images a.draw(renderer) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/matplotlib/artist.py", line 38, in draw_wrapper return draw(artist, renderer, args, kwargs) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/cartopy/mpl/geoaxes.py", line 479, in draw return matplotlib.axes.Axes.draw(self, renderer=renderer, kwargs) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/matplotlib/artist.py", line 38, in draw_wrapper return draw(artist, renderer, args, kwargs) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/matplotlib/axes/_base.py", line 2647, in draw mimage._draw_list_compositing_images(renderer, self, artists) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/matplotlib/image.py", line 135, in _draw_list_compositing_images a.draw(renderer) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/matplotlib/artist.py", line 38, in draw_wrapper return draw(artist, renderer, *args, *kwargs) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/cartopy/mpl/feature_artist.py", line 155, in draw geoms = self._feature.intersecting_geometries(extent) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/cartopy/feature/init.py", line 302, in intersecting_geometries return super(NaturalEarthFeature, self).intersecting_geometries(extent) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/cartopy/feature/init.py", line 110, in intersecting_geometries return (geom for geom in self.geometries() if File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/cartopy/feature/init.py", line 287, in geometries geometries = tuple(shapereader.Reader(path).geometries()) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/cartopy/io/shapereader.py", line 140, in init self._reader = reader = shapefile.Reader(filename) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/shapefile.py", line 553, in init self.load(args[0]) File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/shapefile.py", line 648, in load self.shpHeader() File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/shapefile.py", line 730, in shpHeader self.shpLength = unpack(">i", shp.read(4))[0] 2 struct.error: unpack requires a buffer of 4 bytes

CODE:

import numpy as np import pandas as pd
import json
from os import path
from datetime import datetime, timedelta
import pdb

TAWES_ID_out=[5220108,11315,11121,5220120,5220091,5220092,5221029,5220059,5220107,5220117,5220101,5221021,11112,5220114,5220110,5221005,5220002,8989099,11117,5220051] # INN

pd_TAWES_info = pd.read_csv('TAWES_Sations.csv')

kk = [] # Lat lon data

for ID in TAWES_ID_out:
try:
kk.append(pd_TAWES_info[pd_TAWES_info['StatId'] == ID].values[0])
except:
print('Station not found: %d' %(ID))

lon_min = pd_TAWES_info['StatLon'].min()
lon_max = pd_TAWES_info['StatLon'].max()
lat_min = pd_TAWES_info['StatLat'].min()
lat_max = pd_TAWES_info['StatLat'].max()

lon_min = 8.222777777777777
lon_max = 17.685555555555556
lat_min = 45.55166666666666
lat_max = 49.46916666666666

print('lon_min,lon_max,lat_min,lat_max=,', lon_min,lon_max,lat_min,lat_max)

kk = np.array(kk)
print (kk)

pdb.set_trace()

from matplotlib import pyplot as plt
import matplotlib.colors as colors
from matplotlib.colors import BoundaryNorm
from matplotlib import cm
from pathlib import Path
import cartopy
import cartopy.crs as ccrs
from cartopy.feature import NaturalEarthFeature, ShapelyFeature

import cartopy.io.shapereader as shpreader
import cartopy.feature as cfeature
land_50m = cfeature.NaturalEarthFeature('physical', 'land', '50m',
edgecolor='face',
facecolor=cfeature.COLORS['land'])
resol = '10m'
rivers = cartopy.feature.NaturalEarthFeature( category='physical', name='rivers_lake_centerlines', scale='10m')
border = cartopy.feature.NaturalEarthFeature('cultural', 'admin_0_countries', '10m')
lakes = cartopy.feature.NaturalEarthFeature(category='physical', name='lakes', scale='10m')
lakes = cartopy.feature.NaturalEarthFeature('physical', 'lakes', \
scale=resol, edgecolor='b', facecolor=cfeature.COLORS['water'])
land = cartopy.feature.NaturalEarthFeature('physical', 'land', \
scale=resol, edgecolor='k', facecolor=cfeature.COLORS['land'])
ocean = cartopy.feature.NaturalEarthFeature('physical', 'ocean', \
scale=resol, edgecolor='none', facecolor=cfeature.COLORS['water'])

fig=plt.figure(figsize=(18, 7))
ax = plt.axes(projection=cartopy.crs.Mercator())

ax = plt.subplot(1, 1, 1, projection=cartopy.crs.Mercator())

plt.title('Stations')

MUST BE LON LAT; NOT LAT LON

ax.plot(kk[:,2],kk[:,1],'or',transform=cartopy.crs.PlateCarree())
ax.set_extent([lon_min, lon_max, lat_min, lat_max],crs=cartopy.crs.PlateCarree())

transform = ccrs.PlateCarree()._as_mpl_transform(ax)

for i in range(len(kk)):
ax.annotate('%d' %(kk[i,0]), xy=(kk[i,2],kk[i,1]), xycoords=transform,ha='right', va='top')

ORDER MATTERS;LAND FIRST!

ax.add_feature(land,facecolor='beige')
ax.add_feature(ocean, linewidth=0.2 )
ax.add_feature(border, edgecolor='k', linewidth=2, linestyle='-', facecolor='none')

ax.add_feature(provinces, edgecolor=0.5, linewidth=1, linestyle='-', facecolor='none')

ax.add_feature(rivers, edgecolor='darkblue', linewidth=1, linestyle='-',facecolor='none', alpha=0.5)

ax.add_feature(lakes, edgecolor='darkblue', linewidth=1, linestyle='-',facecolor='none', alpha=0.5)

ax.add_feature(cartopy.feature.LAKES, edgecolor='darkblue', alpha=0.5)

ax.add_feature(lakes)

ax.add_feature(lakes_10m)

ax.add_feature(cfeature.LAKES)

plt.savefig('Station_plot.png')

plt.show()

[greglucas]

This looks like an issue with pyshp (inside the shapefile.py file) not being able to read the file. You can try to install Fiona, which Cartopy will then use to read the shapefiles instead. conda install fiona

[alexandrefierro]

Dear Greg:

I tried the solution you proposed but still no luck (which leads me to this question: why aren't those "shapefiles" simply not written in a more transparent/elegant/workable format such as NETCDF?):

python -m pip install fiona

then got:

File "/mapp_arch3/afierro/anaconda3/envs/env_inca/lib/python3.7/site-packages/fiona/collection.py", line 162, in init self.session.start(self, **kwargs) File "fiona/ogrext.pyx", line 540, in fiona.ogrext.Session.start File "fiona/_shim.pyx", line 90, in fiona._shim.gdal_open_vector fiona.errors.DriverError: '.local/share/cartopy/shapefiles/natural_earth/physical/ne_10m_lakes.shp' not recognized as a supported file format.

[greglucas]

Looks like you have a bad file then. Perhaps try downloading the file manually https://www.naturalearthdata.com/downloads/10m-physical-vectors/ and putting all of the files into the cache location after extracting .local/share/cartopy/shapefiles/natural_earth/physical/.

[alexandrefierro]

After realizing that some of my lake shapefiles at 10 and 50m had issues (for some unknown reasons), I tried to plot the 110m res - the code ran fine this time but no lakes were plotted...: 0 -rw-r--r-- 1 afierro 0 Nov 26 07:39 ne_10m_lakes.shp 4 -rw-r--r-- 1 afierro 330 Nov 29 16:12 ne_10m_land.dbf 6468 -rw-r--r-- 1 afierro 6589552 Nov 29 16:12 ne_10m_land.shp 4 -rw-r--r-- 1 afierro 180 Nov 29 16:12 ne_10m_land.shx 4 -rw-r--r-- 1 afierro 138 Nov 29 16:15 ne_10m_ocean.dbf 7040 -rw-r--r-- 1 afierro 7176016 Nov 29 16:15 ne_10m_ocean.shp 4 -rw-r--r-- 1 afierro 108 Nov 29 16:15 ne_10m_ocean.shx 9616 -rw-r--r-- 1 afierro 9797637 Jun 13 17:35 ne_10m_rivers_lake_centerlines.dbf 4112 -rw-r--r-- 1 afierro 4187664 Jun 13 17:35 ne_10m_rivers_lake_centerlines.shp 12 -rw-r--r-- 1 afierro 11740 Jun 13 17:35 ne_10m_rivers_lake_centerlines.shx 0 -rw-r--r-- 1 afierro 0 Nov 26 10:10 ne_110m_coastline.shp 164 -rw-r--r-- 1 afierro 163315 Nov 29 16:03 ne_110m_lakes.dbf 12 -rw-r--r-- 1 afierro 9324 Nov 29 16:03 ne_110m_lakes.shp 4 -rw-r--r-- 1 afierro 300 Nov 29 16:03 ne_110m_lakes.shx 0 -rw-r--r-- 1 afierro 0 Nov 26 10:10 ne_50m_lakes.shp 76 -rw-r--r-- 1 afierro 69709 Nov 26 10:34 ne_50m_land.dbf 1036 -rw-r--r-- 1 afierro 1050332 Nov 26 10:34 ne_50m_land.shp 12 -rw-r--r-- 1 afierro 11460 Nov 26 10:34 ne_50m_land.shx

[greglucas]

Try setting the zorder of the feature, general questions are better posed on Stackoverflow. Closing this as it looks like this is an issue with your system and not with Cartopy.

[dopplershift]

I tried the solution you proposed but still no luck (which leads me to this question: why aren't those "shapefiles" simply not written in a more transparent/elegant/workable format such as NETCDF?):

We're simply downloading them from the publicly available Natural Earth dataset, which is now hosted on Amazon S3. Feel free to direct any and all comments regarding formats to them . I will say, though, that shapefiles are pretty much an industry standard for storage of vector geographic datasets. As much as I love netCDF, it's generally not particularly well-suited to (or generally used for) this kind of data since individual geographic features have varying numbers of points (as opposed to standard gridded raster data).

[alexandrefierro]

I tried the solution you proposed but still no luck (which leads me to this question: why aren't those "shapefiles" simply not written in a more transparent/elegant/workable format such as NETCDF?):

We're simply downloading them from the publicly available Natural Earth dataset, which is now hosted on Amazon S3. Feel free to direct any and all comments regarding formats to them . I will say, though, that shapefiles are pretty much an industry standard for storage of vector geographic datasets. As much as I love netCDF, it's generally not particularly well-suited to (or generally used for) this kind of data since individual geographic features have varying numbers of points (as opposed to standard gridded raster data).

Howdy, Ryan; Thanks for the detailed input! Interestingly, I have always viewed boundaries or irregularly shaped objects on a map as data that could easily be rasterized [e.g., boolean array] and, hence, stored as gridded structures in an NETCDF file (pixel-wise). Perhaps the issue in doing so would be an increase in file size, which nowadays isn't really a problem; especially if tapes (despite their slow access) are used to store the data. I have always been vouching for one, single, universal binary (NETCDF) and ASCII (CSV) format as juggling between those (e.g., netcdf3,4 <-> grib1,2, grib1<-> grib2, native binaries <-> grib2 etc...) is often non-trivial. As far as computer language go; I would vouch for Fortran across the board because one is inclined to actually build her/his own functions from scratch instead of relying on third-party "black boxes" for computations (e.g., Python libraries). Warm wishes, Alex-

[dopplershift]

Not to get to far in the weeds here but:

1. Converting vector (polygon/line) data to a raster data requires "sampling" the geometry, which introduces some jitter (because you're nearest neighbor sampling to a 0/1 grid). This can be mitigated by having a super fine grid, but that introduces a speed & space vs. accuracy tradeoff. You can likely come up with a good tradeoff for YOUR applications, but good luck coming up with one that suits all, or even most, applications. You're also introducing all sorts of "jaggy" issues that you see in computer graphics and rendering lines and triangles.
2. While storage space is not necessarily an issue, ballooning memory to load that grid and work with it certainly could be. Also, compute speed due to the need to process all the points in a grid rather than only the points defining a geometry is a potential problem, and almost guaranteed to be encountered in many applications.
3. Appropriate gridding for the vector data probably depends on the end projected space you want to work in--just imagine what the appropriate raster grid would be for a polygon near the Earth's poles in various projected coordinate spaces. This is the same as trying to display satellite data in various projections. Raw vector data are much more easily reprojected.
4. For many applications, vector geometry data are much more readily applicable to the problem at hand. Things like "Does this feature (polygon) contain this point?". Reassembling polygon features from a raster grid is a decidedly non-trivial computer problem.

Let's just say there's a reason the thousands of smart people who have developed the GIS software stack (who work with raster data regularly) haven't decided to convert all boundary data to raster grids. 😉 I don't dispute that such a solution could solve your bottlenecks, but it's definitely not a general solution suitable for a general-purpose dataset.