2

Using GeoPandas and Contextily to create maps with a Stamen background map I notice a misalignment between the background map and the coordinates of polygons drawn on the map. The shapefile is from the official Dutch Bureau of Statistics so the coordinates are correct.

What can I do?

My usage should be inline with the instruction and the example notebook. Others seem to have the same problem, e.g. this tutorial.

My original projection is:

{'proj': 'sterea',
 'lat_0': 52.15616055555555,
 'lon_0': 5.38763888888889,
 'k': 0.9999079,
 'x_0': 155000,
 'y_0': 463000,
 'ellps': 'bessel',
 'units': 'm',
 'no_defs': True}

I perform:

dfmap = dfmap.to_crs(epsg=3857)

The crs becomes:

{'init': 'epsg:3857', 'no_defs': True}

I then perform:

ax = df.plot(figsize=(10, 10), alpha=0.5, edgecolor='k')
ctx.add_basemap(ax)

But the basemap and polygon are misaligned, as per the map image (the polygon should be more to the left and a bit lower so that the curve aligns with the water on the map):

example of misaligned map

I tried different zoomlevels, checking axis values, setting ax.axis('square'),...

On axis values, adding the basemap changes the values:

print(ax.axis())
ctx.add_basemap(ax, zoom=15)
print(ax.axis())

Returns:

(542505.2085419805, 543235.6744236776, 6866130.598686902, 6866937.041798945)
(541785.6564853293, 544231.6413904548, 6865879.628687669, 6867102.621140235)

Which I would not expect considering the following in add_basemap(ax)?:

xmin, xmax, ymin, ymax = ax.axis()
# restore original x/y limits
ax.axis((xmin, xmax, ymin, ymax))

I also tried to add the basemap after performing dfmap.to_crs(epsg=4326) but then I get the sea as background. And I also tried with ll=True as well as ll=False in the ctx.bounds2img function:

def add_basemap(ax, zoom, url='http://tile.stamen.com/terrain/tileZ/tileX/tileY.png'):
    xmin, xmax, ymin, ymax = ax.axis()
    basemap, extent = ctx.bounds2img(xmin, ymin, xmax, ymax, zoom=zoom, url=url, ll=True)
    ax.imshow(basemap, extent=extent, interpolation='bilinear')
    # restore original x/y limits
    ax.axis((xmin, xmax, ymin, ymax))

I lack the knowledge about long lat vs xy but if I read about how "the Web Mercator uses the spherical formulas at all scales whereas large-scale Mercator maps normally use the ellipsoidal form of the projection. The discrepancy is imperceptible at the global scale but causes maps of local areas to deviate slightly from true ellipsoidal Mercator maps at the same scale. This deviation becomes more pronounced further from the equator, and can reach as much as 40 km on the ground." then I guess it relates to that.

But how?

  • thanks for the detailed question. I don't directly see something wrong, and normally this should work (Web Mercator should not give such an error at that scale). Would you be able to provide a small sample dataset (to create dfmap to try to reproduce the issue) ? – joris Jan 23 at 12:08
  • Thank you so much! I put the sample information in this online folder: 1drv.ms/u/s!AlNIYisf3zAHs8cBi7cUPP-JS9Yl3Q?e=Iqe5HP – Wouter Jan 23 at 12:40
  • 1
    Looks like a datum transformation issue (or lack thereof). 4326 and 3857 use WGS84 as a GCS. You need to transform the Bessel to that for (closer) alignment. – Vince Jan 23 at 12:50
  • Thanks. So I cannot use the GeoPandas method .to_crs() in this case and need to do the re-projecting from Bessel to EPSG 3857 in an alternative manner? I will investigate how to do that and get back here. – Wouter Jan 24 at 10:22
1

The root cause seems to be an incorrect (or incomplete) definition of the CRS of your data. Your data is actually using EPSG:28992 (https://epsg.io/28992), and if you look at the PROJ.4 definition on that webpage is not fully matching the crs that is returned by the original dfmap.crs (it is missing the "towgs84" part).

When "correcting" that by overwriting the crs attribute with the PROJ.4 string from https://epsg.io/28992, I get the expected result:

dfmap.crs = "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs"
dfmap2 = dfmap.to_crs(epsg=3857)

ax = dfmap2.plot(figsize=(10, 10), alpha=0.5, edgecolor='k')
contextily.add_basemap(ax)
ax.set_axis_off()

enter image description here

Now, as to why this happens for this specific CRS, I don't know. The files you have include a proper WKT description of the CRS in the .prj file of the Shapefile.

  • Also the map from my question is now perfectly aligned, fantastic, many thanks! – Wouter Jan 24 at 14:55
  • 1
    Update: with the development version of geopandas (to become 0.7), it is working correctly out of the box. The main difference is that we are now requiring a more recent version of pyproj, and when reading the .prj information from the file better preserve the information of the CRS (no longer rely on a roundtrip through PROJ4 strings) – joris Jan 25 at 9:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.