(Un)wrapping of DEM tiles across dateline

[cmarshak]

Solves #47 by (un)wrapping DEM tiles and geoid when extents cross anti-meridian / dateline.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[cmarshak]

Here are two examples that include typical resampling for square resolution cells at these higher latitudes:

aleutian_example_l.tif
aleutian_example_r.tif

Unfortunately, QGIS has some difficulty wrapping map tiles at the dateline too. When using epsg:3857 for viewing in QGIS, the rasters wrapped automatically so a lat/lon that extended beyond the dateline appeared as disjoint rasters. When using epsg:4326, the memory on my machine spiked and I had to force quit. In the next comment, I provide some tiles for more direct comparison, albeit much more annoying.

The rasters were generated as indicated here:

from dem_stitcher.stitcher import stitch_dem
import rasterio
from pathlib import Path

# Aleutian Volcano Chain
# xmin, ymin, xmax, ymax
bounds_l = [-181, 51, -179, 52]
bounds_r = [179, 51, 181, 52]

dst_area_or_point = 'Point'
dst_ellipsoidal_height = False
dem_name = 'glo_30'

X_l, p_l = stitch_dem(bounds_l,
                      dem_name=dem_name,
                      dst_ellipsoidal_height=dst_ellipsoidal_height,
                      dst_resolution=0.00027777777777777,
                      dst_area_or_point=dst_area_or_point)

X_r, p_r = stitch_dem(bounds_r,
                      dem_name=dem_name,
                      dst_ellipsoidal_height=dst_ellipsoidal_height,
                      dst_resolution=0.00027777777777777,
                      dst_area_or_point=dst_area_or_point)

with rasterio.open('aleutian_example_l.tif', 'w', **p_l) as ds:
    ds.write(X_l, 1)
    ds.update_tags(AREA_OR_POINT=dst_area_or_point) 

with rasterio.open('aleutian_example_r.tif', 'w', **p_r) as ds:
    ds.write(X_r, 1)
    ds.update_tags(AREA_OR_POINT=dst_area_or_point) 

[cmarshak]

Copying tags from the tiles ensures rest of workflow is correct. If relevant tags are not present, the tiles are assumed to have AREA tag and there will be a nefarious half pixel shift. I realized this when I made an (embarrassing) post on the Affine library and Sean Gillies very nicely responded indicating that I was not writing down the affine transformation correctly.

[cmarshak]

This demonstrates that approaching Aleutians from the left of the globe or the right provides identical arrays and the same geo-metadata after translation through the stitch_dem API.