GeoSensorWebLab / polarmap.js

Custom Leaflet layer for re-projecting maps and map features
BSD 2-Clause "Simplified" License
53 stars 23 forks source link

Track down "magic" numbers #3

Closed openfirmware closed 8 years ago

openfirmware commented 8 years ago

On the LAEA tile configuration, I am using a magic number for the projection's extent:

var extent = 11000000 + 9036842.762 + 667;

I honestly don't know where all of these numbers come from. I added the 667 to get the map to "line up" markers properly as they were offset a bit from their true positions (when compared to the marker positions on an EPSG:3857 map such as OpenStreetMap). 9036842.762 is exactly the extent for the LAEA projections 3571 – 3576. But the 11 million is a mystery.

While the map appears to "work" with the tile backend, I am not sure why these numbers work and I would like to know. If it helps, here is the configuration for the tiles on the server side. The tiles are rendered with Mapnik and renderd/mod_tile.

I may have to set up a local installation of the tile server and PolarMap.js to test changes on both the server and client at the same time. If so, I will write up a guide on how to do that.

timrobertson100 commented 8 years ago

I am also digging into this, and I believe it might relate to the fact that while the recommended envelope for the projection is above 45 degrees, you can paint a tile to world extent with heavy distortion at the edges if you wish.

Therefore you need to consider the world bounds, which means the circumference of the earth at the equator (approx 40 million meters). If you add those numbers up, you get approx 20 million which correlates to 40 million when considering you have +20 and -20 million.

I have written my own data overlays, and it was this logic that I needed to apply in my scaling AffineTransform to align .

openfirmware commented 8 years ago

Interesting! I recall now I had to cut all the shapefiles to exclude data below the equator, as the distortion caused glitchy rendering with Antarctica. Here is what I mean in QGIS:

no-cut-epsg-3573

And with a crop:

cut-epsg-3573

(Scale is 1:120,000,000, EPSG:3573 on-the-fly reprojection)

That would explain the empty space around the world at zoom level 0.

Your suggestion regarding the equator seems to be the solution, and the "9036842.762" number is just me extracting a number that isn't actually relevant.

The LAEA projections of EPSG:3571 to 3576 are all using the WGS84 ellipsoid. The ellipsoid uses an equatorial radius of 6 378 137 metres, which corresponds to a equatorial circumference of 40 075 016.686 metres. Half that is 20 037 508.343 metres. The magic number in the PolarMap.js code is 20 037 509.762 metres: very close!

Your code is similar to how Proj4Leaflet handles the resolutions. I first define the metres per pixel in the tile layer, the use those to build the resolution zoom steps for a custom Map CRS on Leaflet.

Thanks for the help Tim! I will do some testing with replacing the magic number with the half-circumference and see if markers still line up properly. If so I will release a new version of PolarMap with the new number.

timrobertson100 commented 8 years ago

Great. 40,075,016.686 is also the number I used as I found it in Leaflet and some other GIS libs. Good to know that it corresponds with the WGS84 equatorial radius.

openfirmware commented 8 years ago

I adjusted the extent and saw a tiny shift in marker placement:

comparison

This is to be expected and should be a more accurate representation of the location.

duckontheweb commented 8 years ago

I just set these to the extents given by @openfirmware here ([-20037508, 20037508]) without any other adjustments and it seemed to pinpoint the location accurately. I may be misunderstanding things, but I don't think that the only thing that the origin property in proj4leaflet needs to know is where the tile origin point is, and that can come directly from the projected tile extents.

@openfirmware Are the tile extent settings the same for the other Arctic projections? If not, could you point me to the configs for those?

Thanks

openfirmware commented 8 years ago

The extent settings should be the same for the other projections, the only difference in configuration on the tile server (https://github.com/GeoSensorWebLab/awm-styles/blob/master/awm-styles.yaml) is that the Spatial Reference System string has a different central meridian.

duckontheweb commented 8 years ago

Awesome, thanks. After digging into this a bit more, it looks like that origin/extent value comes from the major axis after all. Should have just followed your link in the first place... Thanks for the help.