Match projection of shapefile in cartopy

Question

I am trying to make a Choropleth map using matplotlib and cartopy for which I obviously need to plot a shapefile first. However, I did not manage to do so, even though a similar question has been asked here and here. I suspect either the projection or the bounds to be misspecified.

My shapefile has the projection

PROJCS["WGS_1984_UTM_Zone_32Nz",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",32500000],
    PARAMETER["False_Northing",0],
    PARAMETER["Central_Meridian",9],
    PARAMETER["Scale_Factor",0.9996],
    PARAMETER["Latitude_Of_Origin",0],
    UNIT["Meter",1]]

and can be downloaded here where I am talking about vg250_2010-01-01.utm32w.shape.ebenen/vg250_ebenen-historisch/de1001/vg250_gem.shp

My code is

#!/usr/local/bin/python
# -*- coding: utf8 -*-

import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import matplotlib.pyplot as plt

fname = 'path/vg250_gem.shp'
proj = ccrs.TransverseMercator(central_longitude=0.0,central_latitude=0.0,
                           false_easting=32500000.0,false_northing=0.0,
                           scale_factor=0.9996)
municipalities = list(shpreader.Reader(fname).geometries())
ax = plt.axes(projection=proj)
plt.title('Deutschland')
ax.add_geometries(municipalities,proj,edgecolor='black',facecolor='gray',alpha=0.5)
ax.set_extent([32458044.649189778*0.9, 5556418.748046352*1.1, 32465287.307457082*0.9, 5564153.5456742775*1.1],proj)
plt.show()

where I obtained the bounds using the corresponding method from fiona. Python throws an error

Traceback (most recent call last):
  File "***/src/analysis/test.py", line 16, in <module>
ax.set_extent([32458044.649189778, 5556418.748046352, 32465287.307457082, 5564153.5456742775],proj)
  File "/usr/local/lib/python2.7/site-packages/cartopy/mpl/geoaxes.py", line 652, in set_extent
ylim=self.projection.y_limits))
ValueError: Failed to determine the required bounds in projection coordinates. Check that the values provided are within the valid range (x_limits=[-20000000.0, 20000000.0], y_limits=[-10000000.0, 10000000.0]).
[Finished in 53.9s with exit code 1]

This doesn't make sense to me. Also, experimenting with ccrs.UTM() gives me a plot showing a white area. I'd appreciate it if anyone can tell me how to fix this. Thank you!

I don't know anything about cartopy etc., but reading the next-to-last line in the error message, there are some required x and y limits listed, but the limits in your call to `set_extent` (the call which caused the error) are outside those required limits. — tsj, Oct 08 '16 at 19:32
That is what I referred to with "This doesn't make sense to me" since I extracted the bounds of the shapefile using fiona. Above, I tried experimenting with the `set_extent` method - including deleting this line altogether - but that didn't work either. — Jhonny, Oct 09 '16 at 08:52
Regardless of what bounds you were given, they are outside the bounds that cartopy is willing to accept. Again, I'm not familiar with these packages, but I notice that in the construction of `proj`, you supply a `false_easting` of 32.5M. Presumably this would result in a shift in the x direction. In your call to `set_extent`, you have two values near 32.5M that would be corrected by such a shift. Looking at the docs for cartopy, the input to `set_extent` is (x0,x1,y0,y1), but it looks like you may have given (x0,y0,x1,y1), so the `false_easting` would not get you back to acceptable bounds. — tsj, Oct 09 '16 at 21:54

tsj · Accepted Answer · 2016-10-10T15:51:59.007

I have found two issues. One is an incorrect specification of limits in your call to set_extent, the documentation specifies [x0,x1,y0,y1] should be the input, you seem to have given [x0,y0,x1,y1]. The other issue seems to be a limitation in cartopy, as best I can tell. It looks like projections outside the limits listed in the error message will always fail, and those limits are hardcoded. You can just edit the source (this line in their latest release), changing -2e7 to -4e7, likewise for the upper bound. After these fixes, your plot is generated without issue: The new set_extent line:

ax.set_extent([32458044.649189778*0.975, 32465287.307457082*1.025,5556418.748046352*0.9, 556415,3.5456742775*1.1],proj)

You may also want to set central_longitude=9.0 in your TransverseMercator, that seems to be what's specified in your shapefile.

I would recommend contacting the developers about this, they might have a good reason for setting those bounds, or they might have a better workaround, or maybe they'll widen the bounds in a later release!

Update

Your bounds also seem to have been set based on only the first of the municipalities:

In [34]: municipalities[0].bounds
Out[34]: (32458044.649189778, 5556418.748046352, 32465287.307457082, 5564153.5456742775)

But the other elements have different bounds. You can get limits flushed to the actual drawing based on min/max values of the bounds of all municipalities.

bnd = np.array([i.bounds for i in municipalities])
x0,x1 = np.min(bnd[:,0]),np.max(bnd[:,2])
y0,y1 = np.min(bnd[:,1]),np.max(bnd[:,3])
ax.set_extent([x0,x1,y0,y1],proj)

Thank you for your effort, @tsj. I don't quite understand why the bounds of the shapefile that I got using fiona don't give me exactly the smallest possible square in which the country would fit. Is the projection (still) misspecified in my code from above? Is there a mistake in Fiona, or cartopy? At least an inconsistency? — Jhonny, Oct 10 '16 at 09:35
Also, just a thought, +/- 2e7 should be sufficient from the cartopy side, since the circumference of the earth is about 40e6. Maybe there's some way to change your input so that it's referenced to -7.5e6 instead of 32.5e6, then you wouldn't have to go edit the source of cartopy. — tsj, Oct 10 '16 at 16:01
So the geometry method actually returns a shapley object which is why I can call `bounds`, right? And because this returns `(x0, y0, x1, y1)` while cartopy expects `(x0,x1,y0,y1)` I need to reorganize as you're showing in the code snippet. Got it, thanks - works like a charm. — Jhonny, Oct 13 '16 at 16:23

Match projection of shapefile in cartopy

1 Answers1

Linked