pythongisshapefilemap-projectionsmapnik

Mapnik not -85.0511º to 85.0511º shape issue


I have a shape that boundaries goes from -180º to 180º on longitude, but on latitude, It goes from -90º to 83.64513º. In a front-end application using leaflet, when It asks for mapnik server the image of the tiles, I convert tile positions to latitude/longitude. Longitude works, but latitude doesn't. I'm using this formula to convert:

lat = arctan(sinh(pi*(1 - 2*y/2^zoom))) * 180/pi

The problem (I don't really know if It's the problem...) is that this formula admits that latitude goes from -85.0511º to 85.0511º, then I'm getting erros like this image:

enter image description here

What could I do to workaround this problem? Change the shape size (then how can I do it?), maybe there's a generic formula that I could pass any latitude or I'm missing some step.

Instead of using the formula above, I also tried using this code with the GoogleTile method. Got the same results...

Here is the code I'm using:

@app.route('/tiles/<z>/<x>/<y>', methods=['GET'])
def tiles(z, x, y):
    filename = tiles_path + r"tile_%s_%s_%s.png" % (z, x, y,)
    filename = filename.encode('ascii', 'ignore')
    z = float(z); x = int(x); y = int(y)
    if not os.path.isfile(filename):
        x_1, y_1 = num2deg(x, y, z)
        x_2, y_2 = num2deg(x + 1, y + 1, z)
        envelope = mapnik.Envelope(x_1, y_1, x_2, y_2)
        mapnik_map.zoom_to_box(envelope)
        mapnik.render_to_file(mapnik_map, filename, "png")
    return send_file(filename)

def num2deg(xtile, ytile, zoom):
    n = 2.0 ** zoom
    lon_deg = xtile / n * 360.0 - 180.0
    lat_rad = atan(sinh(pi * (1 - 2 * ytile / n)))
    lat_deg = degrees(lat_rad)
    return lon_deg, lat_deg

If I change the aspect_fix_mode to ADJUST_CANVAS_HEIGHT:

mapnik_map.aspect_fix_mode = mapnik.aspect_fix_mode.ADJUST_CANVAS_HEIGHT

I don't have the problem above, but setting this, I get the map stretched out and distorted.

Here is the shapefile I'm using.


EDIT:

Mapnik default projection:

map_obj.srs
>>> '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
map_obj.layers[0].srs
>>> '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'

EDIT 2:

Tried to modify the map and layer srs like this:

mapnik_map = mapnik.Map(256, 256, "+init=epsg:3857")
...
layer = mapnik.Layer("layer_name", "+init=epsg:4326")

But nothing is rendered specially when the map srs is set.


Solution

  • Well, I was making a huge mistake about which projection standard to use. As my data is in epsg:4326, I decided to change everything to suit this standard. Here is what I did to make things work:

    # creating the map
    map = mapnik.Map(map_size, map_size, '+init=epsg:4326')
    
    # creating a layer
    layer = mapnik.Layer('layer', "+init=epsg:4326")
    
    # tile to degree conversion (globalmaptiles.py adaptation)
    # GlobalGeodetic.TileBounds
    def tile2deg(tx, ty, zoom):
        res = 180 / 256.0 / 2**zoom
        return (
            tx*256*res - 180,
            ty*256*(-res) + 90,
            (tx+1)*256*res - 180,
            (ty+1)*256*(-res) + 90
        )
    
    # tile2deg usage
    map_bounds = tile2deg(x, y, z)
    envelope = mapnik.Envelope(*map_bounds)
    map.zoom_to_box(envelope)
    

    Front-end:

    // leaflet map configuration
    var map = L.map('map', {
        center: [0, 0],
        zoom: 1,
        subdomains: [],
        crs: L.CRS.EPSG4326,
        tms: false,
    });
    

    Hope to help someone newbie like me in the future :)