GeoTIFF / georaster-layer-for-leaflet

Display GeoTIFFs and soon other types of raster on your Leaflet Map
https://geotiff.github.io/georaster-layer-for-leaflet-example/
Apache License 2.0
283 stars 57 forks source link

Bounds incorrect with EPSG:3031 Polar Stereographic #105

Open rockliffelewis opened 1 year ago

rockliffelewis commented 1 year ago

I am working on a leaflet app which is set in EPSG:3031 Antarctic Polar Stereographic. Provided the TIFF is also in 3031 it does display ok, however it does not display the whole image and has some interesting artifacts.

I have worked around the problem by padding the calculated bounds here: https://github.com/GeoTIFF/georaster-layer-for-leaflet/blob/master/src/georaster-layer-for-leaflet.ts#L984

With original bounds: bounds-original

With bounds padded by 0.8: bounds-padded

DanielJDufour commented 1 year ago

Hi, @rockliffelewis . Thanks so much for submitting this. Would you happen to be able to share the tif? Or create a test case in https://github.com/GeoTIFF/georaster-layer-for-leaflet/tree/master/tests ? Supporting rasters displays in polar regions is personally very important to me :-)

rockliffelewis commented 1 year ago

No worries, here is a little test harness that demonstrates the problem: https://jsfiddle.net/3yztdrq8/

rockliffelewis commented 1 year ago

I've had some time to look into the code, and the issue is that you cant really use a two point bounding box to effectively calculate an area in polar stereographic, because a 'box' ([N, W], [S, E]) in this projection is actually a cone.

On top of that, if you take [N,W],[S,E] as (top left, bottom right), well... in this projection that's not always correct because depending where on the map it is it can be (top right, bottom left), or even (bottom right, top left) and leaflet doesn't handle that very well.

I have found it more consistent to use a 4 point polygon, and in this case the bounds of the image do not line up with the latitudes and longitudes anyway and so you need 4 points to actually describe it accurately

DanielJDufour commented 1 year ago

hi @rockliffelewis , thank you for all the info. definitely makes sense. When I wrote the original reprojection code for this library, I definitely missed some things. In theory, reproject-bbox and geo-extent (which is used elsewhere) should be used in the initBounds function.

DanielJDufour commented 1 year ago

Hi, @rockliffelewis , I'm in the progress of improving georaster-layer-for-leaflet's bounds calculation, but I'm running into some weirdness with the origin of the sample COG you referenced in your code. For some reason, I get different origin values depending on what program I am running:

tiffdump

tiffdump ./data/GA4886_VanderfordGlacier_2022_EGM2008_64m-epsg3031.cog | grep 33922
33922 (0x8482) DOUBLE (12) 6<0 0 0 2.40932e+06 -835572 0>

gdalinfo

gdalinfo ./data/GA4886_VanderfordGlacier_2022_EGM2008_64m-epsg3031.cog | grep "Upper Left"
Upper Left  ( 2409289.215, -835539.341) (109d 7'35.40"E, 66d50'16.98"S)

geotiff.js

const origin = image.getOrigin();
[2409321.727264079,-835571.8532756742,0]

Do you know what program was used to create the tiff? Do you know what the origin value should be?

rockliffelewis commented 1 year ago

I reprojected and cogged the original tiff using GDAL, I'll check the original to see if that lines up.

This is the original geotiff: https://public.services.aad.gov.au/datasets/science/voyage_202122020/multibeam/GA-4886/L3/4886_VanderfordGlacier/32Bit_Geotifs/GA4886_VanderfordGlacier_2022_EGM2008_64m.tiff

anton-seaice commented 1 year ago

bremen_sea_ice_conc_2022_9_9.zip

I think I am having the same problem.

When I add this geotiff to my map (it works in QGIS and validates ok with gdal) all it shows is small triangle in a corner of the map. Both the geotiff and my map are using EPSG:3031

image