EPSG:3413 Stereographic North Pole

[noiv]

I'm trying to get daily NASA tiles into Leaflet and got as far as tiles were loaded and zooming worked.

http://noiv.pythonanywhere.com/gibs/

However, the projection doesn't seem to work, at least by using setView(). Trying to isolate the problem I continued at jsfiddle:

http://jsfiddle.net/noiv/NqAcP/

LInking there to recent versions I adapted the proj string from 'stere' to 'sterea', but now the interface doesn't load any tiles at all.

I'm not even sure if this is right approach, would it be better to use the L.Proj.CRS.TMS tilelayers? How would one prepare the parameter arrays for this projection, then?

[perliedman]

Do you have a reference to "daily NASA tiles", since I'm not familiar with them? What projection do they use, etc.

If the tile set uses TMS, you should definitly be using L.Proj.TileLayer.TMS, since Leaflet's builtin TMS support can only handle the spherical mercator projection.

[noiv]

The projection is EPSG:3413 defined by the NSIDC: http://nsidc.org/data/atlas/epsg_3413.html

The server is kind of a mixture of TMS and WMS: https://wiki.earthdata.nasa.gov/display/GIBS/GIBS+Access+Methods

The capabilities http://map1.vis.earthdata.nasa.gov/wmts-arctic/wmts.cgi?SERVICE=WMTS&request=GetCapabilities

The layer is called: MODIS_Terra_CorrectedReflectance_TrueColor

<ows:BoundingBox crs="urn:ogc:def:crs:EPSG::3413">
  <ows:LowerCorner>-4194304 -4194304</ows:LowerCorner>
  <ows:UpperCorner>4194304 4194304</ows:UpperCorner>
</ows:BoundingBox>

Are these the bounds needed? Same order? The first level has 4 tiles, North Pole as center. There are in total 6 zoom levels with these scale denominators:

<0>28571428.5714285746
<1>14285714.2857142873
<2>7142857.1428571437
<3>3571428.5714285718
<4>1785714.2857142859
<5>892857.1428571430

How do I get to resolutions from this?

[noiv]

I rewrote the code http://jsfiddle.net/noiv/NqAcP/14/ using L.Proj.TileLayer.TMS and the proj from spatialreference for proj4js. It seems the bounds are leftX, lowerY, rightX, upperY. Now it loads tiles, but in the wrong y-order. Projection still doesn't work at all.

What else can I do?

crs3413 = new L.Proj.CRS.TMS('EPSG:3413',
  '+proj=sterea +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs',
  [-4194304, -4194304, 4194304, 4194304],
    {resolutions: [
      892857.1428571430, 
      1785714.2857142859, 
      3571428.5714285718, 
      7142857.1428571437, 
      14285714.2857142873, 
      28571428.5714285746
    ]
  }
);

gibsArctic = new L.Proj.TileLayer.TMS('http://map1a.vis.earthdata.nasa.gov/wmts-arctic/wmts.cgi?TIME=2013-09-17&SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=MODIS_Terra_CorrectedReflectance_TrueColor&STYLE=&TILEMATRIXSET=EPSG3413_250m&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}&FORMAT=image%2Fjpeg', 
  crs3413, {
    maxZoom: 6,
    attribution: 'GIBS',
    tileSize: 512,
    noWrap: true,
    continuousWorld: true
});

leafmap = new L.Map('map', {
  crs: crs3413,
  tms: true,
  continuousWorld: true,
  center: [0, 90],
  zoom: 0,
  layers: [gibsArctic]   
});

[perliedman]

It appears that the problem is that you're using TMS. Apparently WMTS (which I don't know much more about than it exists) does not use TMS's row numbering convention, but rather xyz's convention (rows ordered increasing downwards).

Using L.Proj.CRS and L.TileLayer at least gives me a result that looks correct to me: http://jsfiddle.net/NqAcP/15/

[noiv]

Thx, the tiling is correct, but the projection is not good. The North Pole should appear in the middle of the map. I've added a mouseover to check:

http://jsfiddle.net/noiv/NqAcP/

[perliedman]

From some googling, I think WMTS scaleDenominators are not the same as Proj4Leaflet resolution. Proj4Leaflet defines resolution as projection units per pixel (typically meters per pixel), which is the inverse of scale (pixels per projection unit); Leaflet internally uses scale rather than resolution, but it's a matter of preference.

The scaleDenominator in WMTS is apparently calculated from an imaginary DPI value set to 90.7 in the WMTS spec (why you would want to do that is beyond me). See this link for information and formulas that will probably help in calculating scales/resolutions that can be used in Proj4Leaflet: http://gis.stackexchange.com/questions/29671/mathematics-behind-converting-scale-to-resolution

[noiv]

Umm, it look like an earlier post disappeared here. I had found the correct resolution array and can now apply the math even to the 4326 projection with global coverage of same server. There was only a little issue with the tile order, but I could solve that in the source/getTileUrl. Key was to use the definition from spatialreference.com and not the NSIDC's one. But now it seems I can not use two instances of L.Proj.CRS.TMS at the same time, but I'm still checking, might have to do with Leaflet classes or tilelayers.

In case someone else comes along for same server, here the working code:

var crs3413 = new L.Proj.CRS.TMS('EPSG:3413',
     '+proj=sterea +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs',
    [-4194304, -4194304, 4194304, 4194304],
    {resolutions: [8192, 4096, 2048, 1024, 512, 256]}
);

Thanks for your input - sorry for the hazzle.

[perliedman]

A bit unsure of what you mean by "two instances of L.Proj.CRS.TMS at the same time", but Leaflet only supports a single CRS per map - you can't use two at the same time. Adding two L.Proj.TileLayer.TMS instances with different CRS to the same map most certainly will not work, if that is what you mean.