shapefile: 0 features returned

[Marcos Dione]

Warning: I'm kinda cross posting this over here, which is where the shapefiles come from. Also, have in account that I'm completely noob when talking about SRSs and projections; it took me an afternoon just understanding what more or less the SRS in this email are declaring.

I have an issue with rendering osm-carto maps. OSM maps are designed using carto-css, which is in turn converted to Mapnik XML. In my particular case, I then use a python script that in its core does the following calls:

======== >8 ========

m = mapnik.Map(8*256, 8*256)  # metatile size

mapnik.load_map(m, "Elevation.xml", False)  # non-strict

prj = mapnik.Projection(m.srs)

# the following three comments come from the original script I forked from

# get LatLong (EPSG:4326)
l0 = metatile.coords[0]
l1 = metatile.coords[1]

# this is the only time where we convert manually into WebMerc
# Convert to map projection (e.g. mercator co-ords EPSG:900913)
c0 = prj.forward(mapnik.Coord(l0[0], l0[1]))
c1 = prj.forward(mapnik.Coord(l1[0], l1[1]))

bbox = mapnik.Box2d(c0.x, c0.y, c1.x, c1.y)

m.resize(8*256, 8*256)
m.zoom_to_box(bbox)
m.buffer_size = 128

im = mapnik.Image(8*256, 8*256)
mapnik.render(m, im)

======== 8< ========

We'll come back to it in a moment.

osm-carto have recently switched to a set of shapefiles that represent the oceans and seas, which can be found here and here. These shapefiles have the following SRS:

"+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"

The map itself has this one:

<Map srs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over" background-color="#f2efe9">

They're essentially the same, except for "+over". This has been lengthly discussed back when they switched to similar shapefiles back in 2016, but I'm not sure it has an implication in my case.

What I'm seeing is that no data from those shapefiles are rendered. Turning on mapnik's debug logs I see this:

    agg_renderer: Start processing layer=ocean-lz
    agg_renderer: -- datasource=0x1f03ee0
    agg_renderer: -- query_extent=box2d(617305.4404310820391402,5236547.9337608525529504,930391.5082871665945277,5549634.0016169371083379)
    agg_renderer: Start processing style
    shape_index_featureset: 0 features
    agg_renderer: End processing style
    agg_renderer: End layer processing

    agg_renderer: Start processing layer=ocean-lz
    agg_renderer: -- datasource=0x1f03ee0
    agg_renderer: -- query_extent=box2d(460456.6583899010438472,5080310.6479459535330534,1086628.7941020664293319,5706482.7836581189185381)
    agg_renderer: Start processing style
    shape_index_featureset: 0 features
    agg_renderer: End processing style
    agg_renderer: End layer processing

My understanding is that 0 features means "I couldn't find any objects in that datasource within (or touching) the box2d". Now, my first question would be: those coordinates, in what projection are they? what unit?

That's where I come back to the code above. In the discussion from 2016 linked above, the shapefiles' developer said that "you should not attempt to reproject them, not even a simple round trip conversion Mercator->Geographic->Mercator. This will fail since the polygons slightly extend over the 180-degree-meridian." The code seems to be going from one coord system (EPSG:4326) to another (Google/Web Mercator) and back, but I'm so clueless that I can't really tell. Furthermore, my rendering are in the middle of the Mediterranean, and I don't know if that really has an impact in my case or not.

My further questions are:

* Why is this affecting my setup, while the same files and XML works for OSM's servers? BTW, I'm compiling mapnik form source, so it declares to be 4.0. I'll try to test again with 3.12 tomorrow.

* Should I remove those projections in the code?

* Where can I read to better understand all this?

A complete XML minimal example would be:

<?xml version="1.0" encoding="utf-8"?>
<!DOCTYPE Map[]>
<Map srs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over" background-color="#f2efe9">

<Parameters>
  <Parameter name="scale">1</Parameter>
  <Parameter name="metatile">2</Parameter>
  <Parameter name="name"><![CDATA[OpenStreetMap Carto]]></Parameter>
  <Parameter name="description"><![CDATA[A general-purpose OpenStreetMap mapnik style, in CartoCSS]]></Parameter>
  <Parameter name="bounds">-180,-85.05112877980659,180,85.05112877980659</Parameter>
  <Parameter name="center">0,0,4</Parameter>
  <Parameter name="format">png</Parameter>
  <Parameter name="minzoom">0</Parameter>
  <Parameter name="maxzoom">22</Parameter>
</Parameters>

<Style name="ocean-lz" filter-mode="first">
  <Rule>
    <PolygonSymbolizer fill="#aad3df" />
  </Rule>
</Style>
<Layer name="ocean-lz"
  minimum-scale-denominator="750000"
  srs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over">
    <StyleName>ocean-lz</StyleName>
    <Datasource>
       <Parameter name="file"><![CDATA[data/simplified-water-polygons-split-3857/simplified_water_polygons.shp]]></Parameter>
       <Parameter name="type"><![CDATA[shape]]></Parameter>
    </Datasource>
  </Layer>

<Style name="ocean" filter-mode="first">
  <Rule>
    <PolygonSymbolizer fill="#aad3df" />
  </Rule>
</Style>
<Layer name="ocean"
  maximum-scale-denominator="750000"
  srs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over">
    <StyleName>ocean</StyleName>
    <Datasource>
       <Parameter name="file"><![CDATA[data/water-polygons-split-3857/water_polygons.shp]]></Parameter>
       <Parameter name="type"><![CDATA[shape]]></Parameter>
    </Datasource>
  </Layer>

</Map>

I hope this has not been too long.

[Marcos Dione]

Ok, I understand a little more now. The coordinates in mapnik's logs are in meters, and we can use some tools to convert to something we're more used to, like LatLon:

echo 773465.0948891303269193 5392934.5106466235592961 | cs2cs +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +datum=WGS84 +to +proj=latlong +datum=WGS84
6d56'53.359"E   43d31'46.597"N 0.000

This is definitely around the region I want to render.

I wrote a small script to test whether any geometry on the shapefile lies within (actually, intersects) the bounding box, and voilà, there is one:

======== >8 ========

#! /usr/bin/env python3

import fiona
import shapely.geometry

shapefile = fiona.open('data/water-polygons-split-3857/water_polygons.shp')
bbox = shapely.geometry.box(773465.0948891303269193, 5392934.5106466235592961,
                            774076.5911154140485451, 5393546.0068729072809219)
count = 0

for line in shapefile:
    geometry = shapely.geometry.shape(line['geometry'])

    if geometry.intersects(bbox):
        count += 1

print(count)

======== 8< ========

$ ./water_polygons.py
1

I'm more and more baffled...

[Marcos Dione]

I tested over one of the bboxes in the logs I sent, I get 12 hits:

[...]

bbox = shapely.geometry.box( 460456.6583899010438472, 5080310.6479459535330534,
                            1086628.7941020664293319, 5706482.7836581189185381)
[...]

mdione@diablo:~/src/projects/osm$ ./water_polygons.py
12

On Monday, June 3, 2019 at 11:21:40 PM UTC+2, Marcos Dione wrote:

[Marcos Dione]

Even more details: yes, I'm doing a couple of reprojections, but only when calculating the bounding box to rende:

# get LatLong (EPSG:4326)
l0 = metatile.coords[0]
l1 = metatile.coords[1]

# this is the only time where we convert manually into WebMerc
# Convert to map projection (e.g. mercator co-ords EPSG:900913)
c0 = prj.forward(mapnik.Coord(l0[0], l0[1]))
c1 = prj.forward(mapnik.Coord(l1[0], l1[1]))

So I'm going Pixel (in pixels) -> LatLong (in degrees) -> WebMercator (in meters). So that part is fine.

I finally tested mapnik 3.20 again. I renders the shapefiles fine. So the issue applies only to mapnik-master from git. Should I open an issue?

On Monday, June 3, 2019 at 11:21:40 PM UTC+2, Marcos Dione wrote: