[GIS] Using Tilestache and Leaflet To Display Rasters In Native Projection

geotiff-tiffjavascriptleafletmapniktilestache

I have a collection of old Idaho state highway maps in GeoTiff format. The projection that these maps have been georeferenced in is a transverse mercator projection, called IDTM83 (Idaho Transverse Mercator – NAD 83). The projection is defined as:

+proj=tmerc +lat_0=42 +lon_0=-114 +k=0.9996 +x_0=2500000 +y_0=1200000 +ellps=GRS80 +units=m +no_defs

The projection is basically a modified UTM projection that compensates for that fact that the state of Idaho is split by UTM zones 11 and 12. My goal is to serve the maps in the native IDTM83 projection using Tilestache and Leaflet. I am able to serve these same rasters that have already been projected into world mercator, but the the projection distorts much of the line work from the original map.

There have been some helpful resources on the internet on how to serve projected data using Tilestache and Leaflet, here:

  1. https://gist.github.com/skogsmaskin/6452343
  2. http://blog.thematicmapping.org/2012/07/using-custom-projections-with-tilecache.html.

The mapnik configuration for one of my raster maps looks like:

<Map background-color="steelblue" 
srs="+proj=tmerc +lat_0=42 +lon_0=-114 +k=0.9996 +x_0=2500000 +y_0=1200000 +ellps=GRS80 +units=m +no_defs"
maximum-extent="2205517,1134032,2834671,2040323">
<Style name="raster-style">
  <Rule>
      <RasterSymbolizer/>
  </Rule>
</Style>

<Layer name="idaho-1947" status="on" srs="+proj=tmerc +lat_0=42 +lon_0=-114 +k=0.9996 +x_0=2500000 +y_0=1200000 +ellps=GRS80 +units=m +no_defs">
<StyleName>raster-style</StyleName>
<Datasource>
    <Parameter name="type">gdal</Parameter>
    <Parameter name="file">GIS/IdahoMaps/1947_geo_idtm.tif</Parameter>
    <Parameter name="format">tiff</Parameter>
</Datasource>
</Layer>
</Map>

Using Python, I am able to generate maps using bounding boxes from the IDTM coordinate system with the Mapnik API. So I know that Mapnik is not the problem.

My Tilestache config is pretty basic: 

    {
  "cache": {"name": "Test","verbose":"True"},
  "layers": {
    "idaho-1947": {
        "provider": {"name": "mapnik", "mapfile": "./mapnik.xml"},
        "projection": "TileStache.Goodies.Proj4Projection:Proj4Projection('+proj=tmerc +lat_0=42 +lon_0=-114 +k=0.9996 +x_0=2500000 +y_0=1200000 +ellps=GRS80 +units=m +no_defs',[5545984, 2772992, 1386496, 693248, 346624, 173312, 86656, 43328, 21664, 10832, 5416, 2708, 1354, 677, 338.5, 169.25, 84.625, 42.3125, 21.15625, 10.578125, 5.2890625, 1])"
    }}

My leaflet viewer is very basic, set up like on Link #2 above:

  var crs = L.CRS.proj4js('SR-ORG:6677', '+proj=tmerc +lat_0=42 +lon_0=-114 +k=0.9996 +x_0=2500000 +y_0=1200000 +ellps=GRS80 +units=m +no_defs', new L.Transformation(1, -432000, -1, 2040323));
  crs.scale = function(zoom) {
      return 1 / (234.375 / Math.pow(2, zoom));
  };

  var map = new L.Map('map', {
    crs: crs,
    layers: [
      new L.TileLayer('http://localhost:8080/idaho-1947/{z}/{x}/{y}.png', {
        minZoom: 0,
        maxZoom: 4,
        continuousWorld: true
      })
    ],
    center: new L.LatLng(45, -116),
    zoom: 3,
    continuousWorld: true
  });

I intend to refine this once I can get the maps to display, but the example is using UTM data and the scales should be pretty similar. I am stuck here. Using the JavaScript examples provided in both the links above, I cannot display any of the maps – only the lightblue background. I know what I am trying to do is possible in ArcGIS Server, but I would really like to leverage some of the OpenSource tools.

Best Answer

In Tilestache you specify the following array of resolution in the projection object: 

[5545984, 2772992, 1386496, 693248, 346624, 173312, 86656, 43328, 21664, 10832, 5416, 2708, 1354, 677, 338.5, 169.25, 84.625, 42.3125, 21.15625, 10.578125, 5.2890625, 1]

You told us that you are using the same resolutions as link #2 which are:

[234.375, 117.1875, 58.59375, 29.296875, 14.6484375]

This means that your leaflet app probably doesn't query the right tiles.

I looked a little bit at the Leaflet code from your comment. It seems that you still have a different scale array from your TileStache config. A wild guess (untested) would be:

var resolutions = [
  5545984, 2772992, 1386496, 693248, 346624,
  173312, 86656, 43328, 21664, 10832, 5416,
  2708, 1354, 677, 338.5, 169.25, 84.625,
  42.3125, 21.15625, 10.578125, 5.2890625, 1];

// I'm really not sure about this Transformation clause in your case...
var crs = L.CRS.proj4js('SR-ORG:6677', '+proj=tmerc +lat_0=42 +lon_0=-114 +k=0.9996 +x_0=2500000 +y_0=1200000 +ellps=GRS80 +units=m +no_defs'/*, new L.Transformation(1, -432000, -1, 2040323)*/);
crs.scale = function(zoom) {
    return 1 / resolutions[zoom];
};

var map = new L.Map('map', {
  crs: crs,
  layers: [
    new L.TileLayer('http://localhost:8080/idaho-1947/{z}/{x}/{y}.png', {
      minZoom: 0,
      maxZoom: resolutions.length-1,
      continuousWorld: true
    })
  ],
  center: new L.LatLng(45, -116),
  zoom: 3,
  continuousWorld: true
});

If this still doesn't work, you will have to check what are the tile requested. In Chrome, open the dev tool (F12) and look at the Network tab to see what are the tiles requested.

Related Question