Creating an Interrupted Map With OGR, Mapnik and PIL

Question

I am trying to create a simple interrupted Sinusoidal map with Python using OGR, Mapnik and PIL.

The general idea is:

Make 6x 60 degree sinusoidal maps
For each one, make a mask to make it gore shaped (the mask is written to Shapefile for Mapnik)
Glue them together using PIL

So far, I can make a single gore, centred around 0 degrees longitude, but I cannot translate this into the other 6 gores.

~~Gore~~ ^{_{Dead Image Link}}

Obviously, simply changing the area mapped does not give the desired effect:

~~Gore~~ ^{_{Dead Image Link}}

I think that I need to move the central meridian in the sinusoidal projection definition, but when I do so, it just screws up the map:

~~Gore~~ ^{_{Dead Image Link}}

My problem is therefore that I cannot seem to make the other gores. I have tried setting the central meridian as is seen in the Python script, and manually by editing the Proj.4 string in the mapnik XML file, but the above problem remains.

On a side note - if I print out the proj.4 string after importing my edited wkt (as seen in the code below), it always comes out as:

+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

It looks like this means that OGR is not interpreting the central meridian value correctly?

I've searched quite extensively, and have been unable to figure out what the problem is - does anyone know how I can do this?

Here is my Python code (currently only makes one gore, not 6):

#!/usr/bin/env python

import mapnik, ogr, osr, pyproj, os
from PIL import Image

###
# Draw a Rhumb line with nPoints nodes
###
def getRhumb(startlong, startlat, endlong, endlat, nPoints):

    # calculate distance between points
    g = pyproj.Geod(ellps='WGS84')

    # calculate line string along path with segments <= 1 km
    lonlats = g.npts(startlong, startlat, endlong, endlat, nPoints)

    # npts doesn't include start/end points, so prepend/append them and return
    lonlats.insert(0, (startlong, startlat))
    lonlats.append((endlong, endlat))
    return lonlats

###
# Write a geometry to a Shapefile
###
def makeShapefile(geom): 

    # set up the shapefile driver
    driver = ogr.GetDriverByName("ESRI Shapefile")

    # remove old shapefile if required     
    if os.path.exists("tmp_mask.shp"):
        driver.DeleteDataSource("tmp_mask.shp")

    # create the data source
    data_source = driver.CreateDataSource("tmp_mask.shp")

    # create the spatial reference, WGS84
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)

    # create the layer
    layer = data_source.CreateLayer("mask", srs, ogr.wkbPolygon)

    # create the feature
    feature = ogr.Feature(layer.GetLayerDefn())

    # Set the feature geometry using the point
    feature.SetGeometry(geom)

    # Create the feature in the layer (shapefile)
    layer.CreateFeature(feature)

    # Destroy the feature to free resources
    feature.Destroy()

    # Destroy the data source to free resources
    data_source.Destroy()


### CODE FROM HERE ###


# USER SETTINGS
CENTRE_MERIDIAN = 0
GORE_WIDTH = 60 # must multiply into 360

# WGS84
source = osr.SpatialReference()
source.ImportFromEPSG(4326)

# World Sinusoidal
wkt_proj = 'PROJCS["JJH Sinusoidal", GEOGCS["GCS_WGS_1984", DATUM["WGS_1984", SPHEROID["WGS_1984",6378137,298.257223563]], PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]], PROJECTION["Sinusoidal"], PARAMETER["False_Easting",0], PARAMETER["False_Northing",0], PARAMETER["Central_Meridian",' + \
str(CENTRE_MERIDIAN) + '], UNIT["Meter",1]]'
target = osr.SpatialReference()
# target.ImportFromEPSG(54008) 
target.ImportFromWkt(wkt_proj)
print target.ExportToProj4() # this prints "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" so central meridian not being set?

# get the transformation to go from WGS84 to Sinusoidal
transform = osr.CoordinateTransformation(source, target)

# this is used to calculate the mainPoints
halfWidth = GORE_WIDTH / 2

# get the main points of the area of interest and transform
mainPoints = ogr.Geometry(ogr.wkbLinearRing)
# mainPoints.AddPoint(CENTRE_MERIDIAN, 90)
# mainPoints.AddPoint(CENTRE_MERIDIAN - halfWidth, 0)
# mainPoints.AddPoint(CENTRE_MERIDIAN, -90)
# mainPoints.AddPoint(CENTRE_MERIDIAN + halfWidth, 0)
mainPoints.AddPoint(0, 90)
mainPoints.AddPoint(-halfWidth, 0)
mainPoints.AddPoint(0, -90)
mainPoints.AddPoint(halfWidth, 0)

# make the gore as 2 lines then combine (using mainPoints in their wgs84 form)
gore = getRhumb(mainPoints.GetX(1), mainPoints.GetY(0), mainPoints.GetX(1), mainPoints.GetY(2), 100) # get the first rhumb (N-S)     
gore2 = getRhumb(mainPoints.GetX(3), mainPoints.GetY(2), mainPoints.GetX(3), mainPoints.GetY(0), 100) # get the second rhumb (S-N)
gore.extend(gore2)

# transform the mainPoints now (they will be Sinusoidal from here)
mainPoints.Transform(transform)

# create ring for the gore
ring = ogr.Geometry(ogr.wkbLinearRing)
for p in gore:
    ring.AddPoint(p[0], p[1])

# create polygon for the mask
clipper = ogr.Geometry(ogr.wkbPolygon)
clipper.AddGeometry(ring)
clipper.CloseRings()

# build mask by differencing the gore from the wkt bounding box and write to shapefile
wkt1 = "POLYGON ((180 -90,176.04395604395603 -90.0,172.087912087912088 -90.0,168.131868131868146 -90.0,164.175824175824175 -90.0,160.219780219780205 -90.0,156.263736263736263 -90.0,152.307692307692321 -90.0,148.35164835164835 -90.0,144.39560439560438 -90.0,140.439560439560438 -90.0,136.483516483516496 -90.0,132.527472527472526 -90.0,128.571428571428555 -90.0,124.615384615384613 -90.0,120.659340659340657 -90.0,116.703296703296701 -90.0,112.747252747252745 -90.0,108.791208791208788 -90.0,104.835164835164832 -90.0,100.879120879120876 -90.0,96.92307692307692 -90.0,92.967032967032964 -90.0,89.010989010989007 -90.0,85.054945054945051 -90.0,81.098901098901095 -90.0,77.142857142857139 -90.0,73.186813186813183 -90.0,69.230769230769226 -90.0,65.27472527472527 -90.0,61.318681318681314 -90.0,57.362637362637358 -90.0,53.406593406593402 -90.0,49.45054945054946 -90.0,45.494505494505489 -90.0,41.538461538461547 -90.0,37.582417582417577 -90.0,33.626373626373635 -90.0,29.670329670329664 -90.0,25.714285714285722 -90.0,21.758241758241752 -90.0,17.80219780219781 -90.0,13.84615384615384 -90.0,9.890109890109898 -90.0,5.934065934065927 -90.0,1.978021978021985 -90.0,-1.978021978021985 -90.0,-5.934065934065927 -90.0,-9.890109890109898 -90.0,-13.84615384615384 -90.0,-17.80219780219781 -90.0,-21.758241758241752 -90.0,-25.714285714285722 -90.0,-29.670329670329664 -90.0,-33.626373626373635 -90.0,-37.582417582417577 -90.0,-41.538461538461547 -90.0,-45.494505494505489 -90.0,-49.45054945054946 -90.0,-53.406593406593402 -90.0,-57.362637362637372 -90.0,-61.318681318681314 -90.0,-65.274725274725284 -90.0,-69.230769230769226 -90.0,-73.186813186813197 -90.0,-77.142857142857167 -90.0,-81.098901098901081 -90.0,-85.054945054945051 -90.0,-89.010989010989022 -90.0,-92.967032967032992 -90.0,-96.923076923076906 -90.0,-100.879120879120876 -90.0,-104.835164835164846 -90.0,-108.791208791208817 -90.0,-112.74725274725273 -90.0,-116.703296703296701 -90.0,-120.659340659340671 -90.0,-124.615384615384642 -90.0,-128.571428571428555 -90.0,-132.527472527472526 -90.0,-136.483516483516496 -90.0,-140.439560439560466 -90.0,-144.39560439560438 -90.0,-148.35164835164835 -90.0,-152.307692307692321 -90.0,-156.263736263736291 -90.0,-160.219780219780205 -90.0,-164.175824175824175 -90.0,-168.131868131868146 -90.0,-172.087912087912116 -90.0,-176.04395604395603 -90.0,-180 -90,-180.0 -86.086956521739125,-180.0 -82.173913043478265,-180.0 -78.260869565217391,-180.0 -74.347826086956516,-180.0 -70.434782608695656,-180.0 -66.521739130434781,-180.0 -62.608695652173914,-180.0 -58.695652173913047,-180.0 -54.782608695652172,-180.0 -50.869565217391305,-180.0 -46.956521739130437,-180.0 -43.043478260869563,-180.0 -39.130434782608695,-180.0 -35.217391304347828,-180.0 -31.304347826086953,-180.0 -27.391304347826086,-180.0 -23.478260869565219,-180.0 -19.565217391304344,-180.0 -15.652173913043484,-180.0 -11.739130434782609,-180.0 -7.826086956521735,-180.0 -3.913043478260875,-180 0,-180.0 3.913043478260875,-180.0 7.826086956521735,-180.0 11.739130434782609,-180.0 15.652173913043484,-180.0 19.565217391304344,-180.0 23.478260869565219,-180.0 27.391304347826093,-180.0 31.304347826086953,-180.0 35.217391304347828,-180.0 39.130434782608688,-180.0 43.043478260869563,-180.0 46.956521739130437,-180.0 50.869565217391312,-180.0 54.782608695652186,-180.0 58.695652173913032,-180.0 62.608695652173907,-180.0 66.521739130434781,-180.0 70.434782608695656,-180.0 74.34782608695653,-180.0 78.260869565217405,-180.0 82.173913043478251,-180.0 86.086956521739125,-180 90,-176.04395604395603 90.0,-172.087912087912088 90.0,-168.131868131868146 90.0,-164.175824175824175 90.0,-160.219780219780205 90.0,-156.263736263736263 90.0,-152.307692307692321 90.0,-148.35164835164835 90.0,-144.39560439560438 90.0,-140.439560439560438 90.0,-136.483516483516496 90.0,-132.527472527472526 90.0,-128.571428571428555 90.0,-124.615384615384613 90.0,-120.659340659340657 90.0,-116.703296703296701 90.0,-112.747252747252745 90.0,-108.791208791208788 90.0,-104.835164835164832 90.0,-100.879120879120876 90.0,-96.92307692307692 90.0,-92.967032967032964 90.0,-89.010989010989007 90.0,-85.054945054945051 90.0,-81.098901098901095 90.0,-77.142857142857139 90.0,-73.186813186813183 90.0,-69.230769230769226 90.0,-65.27472527472527 90.0,-61.318681318681314 90.0,-57.362637362637358 90.0,-53.406593406593402 90.0,-49.45054945054946 90.0,-45.494505494505489 90.0,-41.538461538461547 90.0,-37.582417582417577 90.0,-33.626373626373635 90.0,-29.670329670329664 90.0,-25.714285714285722 90.0,-21.758241758241752 90.0,-17.80219780219781 90.0,-13.84615384615384 90.0,-9.890109890109898 90.0,-5.934065934065927 90.0,-1.978021978021985 90.0,1.978021978021985 90.0,5.934065934065927 90.0,9.890109890109898 90.0,13.84615384615384 90.0,17.80219780219781 90.0,21.758241758241752 90.0,25.714285714285722 90.0,29.670329670329664 90.0,33.626373626373635 90.0,37.582417582417577 90.0,41.538461538461547 90.0,45.494505494505489 90.0,49.45054945054946 90.0,53.406593406593402 90.0,57.362637362637372 90.0,61.318681318681314 90.0,65.274725274725284 90.0,69.230769230769226 90.0,73.186813186813197 90.0,77.142857142857167 90.0,81.098901098901081 90.0,85.054945054945051 90.0,89.010989010989022 90.0,92.967032967032992 90.0,96.923076923076906 90.0,100.879120879120876 90.0,104.835164835164846 90.0,108.791208791208817 90.0,112.74725274725273 90.0,116.703296703296701 90.0,120.659340659340671 90.0,124.615384615384642 90.0,128.571428571428555 90.0,132.527472527472526 90.0,136.483516483516496 90.0,140.439560439560466 90.0,144.39560439560438 90.0,148.35164835164835 90.0,152.307692307692321 90.0,156.263736263736291 90.0,160.219780219780205 90.0,164.175824175824175 90.0,168.131868131868146 90.0,172.087912087912116 90.0,176.04395604395603 90.0,180 90,180.0 86.086956521739125,180.0 82.173913043478265,180.0 78.260869565217391,180.0 74.347826086956516,180.0 70.434782608695656,180.0 66.521739130434781,180.0 62.608695652173914,180.0 58.695652173913047,180.0 54.782608695652172,180.0 50.869565217391305,180.0 46.956521739130437,180.0 43.043478260869563,180.0 39.130434782608695,180.0 35.217391304347828,180.0 31.304347826086953,180.0 27.391304347826086,180.0 23.478260869565219,180.0 19.565217391304344,180.0 15.652173913043484,180.0 11.739130434782609,180.0 7.826086956521735,180.0 3.913043478260875,180 0,180.0 -3.913043478260875,180.0 -7.826086956521735,180.0 -11.739130434782609,180.0 -15.652173913043484,180.0 -19.565217391304344,180.0 -23.478260869565219,180.0 -27.391304347826093,180.0 -31.304347826086953,180.0 -35.217391304347828,180.0 -39.130434782608688,180.0 -43.043478260869563,180.0 -46.956521739130437,180.0 -50.869565217391312,180.0 -54.782608695652186,180.0 -58.695652173913032,180.0 -62.608695652173907,180.0 -66.521739130434781,180.0 -70.434782608695656,180.0 -74.34782608695653,180.0 -78.260869565217405,180.0 -82.173913043478251,180.0 -86.086956521739125,180 -90))"

poly1 = ogr.CreateGeometryFromWkt(wkt1)
difference = poly1.Difference(clipper)
makeShapefile(difference)

# clean up
clipper.Destroy()
difference.Destroy()

# make the map
map = mapnik.Map(500, 500)
mapnik.load_map(map, "Sinusoidal.xml")

# this grows the image if the map dimensions do not fit the canvas dimensions
map.aspect_fix_mode = mapnik.aspect_fix_mode.GROW_CANVAS
map.zoom_to_box(mapnik.Envelope(mainPoints.GetX(1), mainPoints.GetY(0), mainPoints.GetX(3), mainPoints.GetY(2)))

# render to file (and show me it)
mapnik.render_to_file(map, "gores.png")
im = Image.open('gores.png')  
im.show()

(excuse the large wkt statement in the middle, it is used to make the mask).

Here is the mapnik stylesheet:

<?xml version="1.0" encoding="utf-8"?>
<!DOCTYPE Map[]>
<Map srs="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" background-color="#ffffff" >

  <Style name="ne110mland" filter-mode="first">
  <Rule>
    <PolygonSymbolizer fill="#000000" />
  </Rule>
</Style>
<Layer name="ne110mland" srs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs">
    <StyleName>ne110mland</StyleName>
    <Datasource>
       <Parameter name="file"><![CDATA[/.../ne_110m_land/ne_110m_land.shp]]></Parameter>
       <Parameter name="type"><![CDATA[shape]]></Parameter>
    </Datasource>
  </Layer>

  <Style name="mask" filter-mode="first">
  <Rule>
    <PolygonSymbolizer fill="#ffffff" />
    <LineSymbolizer stroke="#000000" stroke-width="4" />
  </Rule>
</Style>
<Layer name="mask" srs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs">
    <StyleName>mask</StyleName>
    <Datasource>
       <Parameter name="file"><![CDATA[/.../tmp_mask.shp]]></Parameter>
       <Parameter name="type"><![CDATA[shape]]></Parameter>
    </Datasource>
  </Layer>

</Map>

jonnyhuck · Accepted Answer · 2016-06-28 17:53:55Z

I have solved it - the map 'screwing up' was, of course, the land masses going beyond +/- 180 degrees. The solution is therefore to clip the land mass data against the gore shape before rendering, which prevents this from happening (and is a little neater than the mask anyway).

I never found out why the projection wasn't updating after I had edited the WKT, but the same problem doesn't seem to occur if you use a Proj.4 string instead.

Corrected code below:

#!/usr/bin/env python

import mapnik, ogr, osr, pyproj, os
from PIL import Image


###
# Draw a Rhumb line with nPoints nodes
###
def getRhumb(startlong, startlat, endlong, endlat, nPoints):

    # calculate distance between points
    g = pyproj.Geod(ellps='WGS84')

    # calculate line string along path with segments <= 1 km
    lonlats = g.npts(startlong, startlat, endlong, endlat, nPoints)

    # npts doesn't include start/end points, so prepend/append them and return
    lonlats.insert(0, (startlong, startlat))
    lonlats.append((endlong, endlat))
    return lonlats

###
# Write a geometry to a Shapefile
###
def makeShapefile(geom, name, layer_name): 

    # set up the shapefile driver
    driver = ogr.GetDriverByName("ESRI Shapefile")

    # remove old shapefile if required     
    if os.path.exists(name):
        driver.DeleteDataSource(name)

    # create the data source
    data_source = driver.CreateDataSource(name)

    # create the spatial reference, WGS84
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)

    # create the layer
    layer = data_source.CreateLayer(layer_name, srs, ogr.wkbPolygon)

    # create the feature
    feature = ogr.Feature(layer.GetLayerDefn())

    # Set the feature geometry using the point
    feature.SetGeometry(geom)

    # Create the feature in the layer (shapefile)
    layer.CreateFeature(feature)

    # Destroy the feature to free resources
    feature.Destroy()

    # Destroy the data source to free resources
    data_source.Destroy()


# USER SETTINGS
CENTRE_MERIDIAN = -150
GORE_WIDTH = 60 # must multiply into 360

# WGS84
source = osr.SpatialReference()
source.ImportFromEPSG(4326)

# World Sinusoidal
target = osr.SpatialReference()
target.ImportFromProj4('+proj=sinu +lon_0=' + str(CENTRE_MERIDIAN) + ' +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs') 

# this is used to calculate the mainPoints
halfWidth = GORE_WIDTH / 2

# get the main points of the area of interest and transform
mainPoints = ogr.Geometry(ogr.wkbLinearRing)
mainPoints.AddPoint(CENTRE_MERIDIAN, 90)
mainPoints.AddPoint(CENTRE_MERIDIAN - halfWidth, 0)
mainPoints.AddPoint(CENTRE_MERIDIAN, -90)
mainPoints.AddPoint(CENTRE_MERIDIAN + halfWidth, 0)

# make the gore as 2 lines then combine (using mainPoints in their wgs84 form)
gore = getRhumb(mainPoints.GetX(1), mainPoints.GetY(0), mainPoints.GetX(1), mainPoints.GetY(2), 100) # get the first rhumb (N-S)     
gore2 = getRhumb(mainPoints.GetX(3), mainPoints.GetY(2), mainPoints.GetX(3), mainPoints.GetY(0), 100) # get the second rhumb (S-N)
gore.extend(gore2)

# create ring for the gore
ring = ogr.Geometry(ogr.wkbLinearRing)
for p in gore:
    ring.AddPoint(p[0], p[1])

# create polygon for the gore
clipper = ogr.Geometry(ogr.wkbPolygon)
clipper.AddGeometry(ring)
clipper.CloseRings()

# write to shapefile
makeShapefile(clipper, "tmp_gore.shp", "gore")

# open countries file and get all of the geometry
shapefile = "ne_110m_land/ne_110m_land.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()
land = ogr.Geometry(ogr.wkbGeometryCollection)
for feature in layer:
    land.AddGeometry(feature.GetGeometryRef())

# clip against the gore
landPanel = land.Intersection(clipper)

# write to shapefile
makeShapefile(landPanel, "tmp_land.shp", "land")

# clean up
clipper.Destroy()
landPanel.Destroy()

# transform the mainPoints now (they will be Sinusoidal from here)
transform = osr.CoordinateTransformation(source, target)
mainPoints.Transform(transform)

# make the map
map = mapnik.Map(500, 500)
mapnik.load_map(map, "Sinusoidal.xml")
map.srs = target.ExportToProj4()

# this grows the image if the map dimensions do not fit the canvas dimensions
map.aspect_fix_mode = mapnik.aspect_fix_mode.GROW_CANVAS

# Set the extent (need to set this to around 0 post transformation as this is the central meridian)
# map.zoom_to_box(mapnik.Envelope(mainPoints.GetX(1), mainPoints.GetY(0), mainPoints.GetX(3), mainPoints.GetY(2)))
map.zoom_to_box(mapnik.Envelope(-3339584.7238, 10001965.7293, 3339584.7238, -10001965.7293))

# render to file (and show me it)
mapnik.render_to_file(map, "gores.png")
im = Image.open('gore.png')  
im.show()

And the associated Manik XML file:

<?xml version="1.0" encoding="utf-8"?>
<!DOCTYPE Map[]>
<Map background-color="#ffffff" >

<Style name="gore" filter-mode="first">
  <Rule>
    <LineSymbolizer stroke="#000000" stroke-width="4" />
  </Rule>
</Style>
<Layer name="gore" srs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs">
    <StyleName>gore</StyleName>
    <Datasource>
       <Parameter name="file"><![CDATA[/.../tmp_gore.shp]]></Parameter>
       <Parameter name="type"><![CDATA[shape]]></Parameter>
    </Datasource>
  </Layer>

  <Style name="ne110mland" filter-mode="first">
  <Rule>
    <PolygonSymbolizer fill="#000000" />
  </Rule>
</Style>
<Layer name="ne110mland" srs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs">
    <StyleName>ne110mland</StyleName>
    <Datasource>
       <Parameter name="file"><![CDATA[/.../tmp_land.shp]]></Parameter>
       <Parameter name="type"><![CDATA[shape]]></Parameter>
    </Datasource>
  </Layer>

</Map>

The completed code is available at github.com/jonnyhuck/GlobeMaker if it is of any use to anyone. — jonnyhuck, Commented Jun 30, 2016 at 12:11

Stack Exchange Network

Creating an Interrupted Map With OGR, Mapnik and PIL

1 Answer 1

Your Answer

Not the answer you're looking for? Browse other questions tagged
python
coordinate-system
ogr
mapnik
or ask your own question.

Hot Network Questions

Creating an Interrupted Map With OGR, Mapnik and PIL

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

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

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
python
coordinate-system
ogr
mapnik
or ask your own question.