Rasterio: Merging rasters that cross the Anti Meridian fails (Goes the long way around)

[ciaranevans]

Expected behavior and actual behavior.

Expected:
Merging 3 rasters results in a merged raster.
Actual:

Unable to allocate array with shape (1, 3875, 1291898) and data type int16: MemoryError
Traceback (most recent call last):
  File "/var/task/lambda_function.py", line 169, in handler
    demall, afdemall = merge.merge(dem_cropped_datasets)
  File "/opt/python/rasterio/merge.py", line 141, in merge
    dest = np.zeros((output_count, output_height, output_width), dtype=dtype)
MemoryError: Unable to allocate array with shape (1, 3875, 1291898) and data type int16

It appears that Rasterio's merge.merge() method doesn't work well with merging tiles which cross over the International Dateline. My guess is because the left most longitude of the tiles is around +178 and the right most longitude of the tiles is around -179. Which is understandably not a usual situation as you'd normally expect L < R.

It looks like Rasterio is trying to create a merged dataset which wraps the longest way around the earth, rather than the IDL.

I assume the code for merge could be altered to account for the IDL and take the shortest route when creating the array for merging?

Steps to reproduce the problem.

dataset_left_of_idl = rasterio.open('path/to/dataset_left.tif', 'r')
dataset_right_of_idl = rasterio.open('path/to/dataset_right.tif', 'r')
datasets_to_merge = [dataset_left_of_idl, dataset_right_of_idl]
merged_datasets, merged_datasets_affine = merge.merge(datasets_to_merge)

I'm using SRTM Tiles from https://e4ftl01.cr.usgs.gov/, you need an account to access it, but if you do, the specific tiles I'm using are:

https://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11/S17E179.SRTMGL1.hgt.zip
https://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11/S17W180.SRTMGL1.hgt.zip
https://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11/S16W180.SRTMGL1.hgt.zip

They cover the North East of Fiji.

Operating system

AWS Lambda (With 3GB max memory)

Rasterio version and provenance

Rasterio 1.1.0 via pip
GDAL 2.4.2

[ciaranevans]

@sgillies many thanks for the label 😄

What is the normal time period for fixes potentially going in?

I'd be interested in contributing a fix but unsure where to start within the merge function... If you've got a rough idea I'd be up for taking a crack at this!

[sgillies]

@ciaranevans here's the thing, which may be a blocker: rasterio's merge requires all inputs to have the same coordinate reference system and requires the output to have that same coordinate reference system. Unless your inputs are in a coordinate reference system centered on the antimeridian, we can't make an output centered on the antimeridian.

A workaround/hack that you could try: wrap the georeferencing of your western hemisphere inputs around so all your files have positive longitudes and merge them to get an output centered on 180 degrees east.

[ciaranevans]

Woops wrong button.

@sgillies thanks for the reply! Ah right, I thought it might be a bit of a PITA!

I took a chance and calculated the distance between the left most bound of the left tile and the right most bound of the right tile and then set the bounds I passed into merge as
left, bottom, left + distance_between_left_and_right, top

It passed and wrote a tiff but the data inside didn't look right. Almost looked empty, I'll give it another go just to make sure I'm not doing something silly there. Otherwise I'll see if I can re-project them all into a joint CRS and see if that works!

It's unfortunate that in the dataset we're trying to generate, there's only that one place in Fiji that falls over haha!

[ciaranevans]

@sgillies tweet was made before I raised this, removed, will be giving your workaround a go!

[ciaranevans]

So, with a bit more Googling and experimenting, we've found a work around @sgillies

Steps go as follows:

1. For each tile to merge calculate_default_transform() to required CRS (In this case, we're going from the SRTM tiles CRS to the MGRS squares
2. Reproject the tile using the calculated values, but pass in the following kwargs:

{'CENTER_LONG': 180}

Appears that https://trac.osgeo.org/gdal/wiki/ConfigOptions#CENTER_LONG will alter the origin of your longitudes

1. Use rasterio.merge() as normal

Gives me:

I think the gaps are due to the reprojection. There is no tile available for the top left so that gap is expected.

[sgillies]

Fiji is such a good test case. And I'm glad to know about CENTER_LONG, that's an obscure one.

There's a lot of overlap beween warping and merging; the latter is kind of a short cut when the input and output have the same CRS. I'm curious about what would happen if you skipped MGRS and called reproject using the same src_crs and dst_crs but with CENTER_LONG set to 180.0, with a destination file that had its left bounds at 179 degrees (say) and right bounds at 181 degrees. Would you be willing to try that?

[ciaranevans]

Yeah, annoyingly the GDAL docs don't elaborate on what it actually is 😂

Just so I'm sure I understand your question, you'd like me to try the reproject with 180 center using the same CRS' for src and dst and then merge after? Happy to give it a go!

[sgillies]

Or even skip the merge. If you create an output dataset that spans the inputs and has a central longitude of 180, you can reproject the inputs, one-by-one, directly into it*.

* Maybe, I haven't tried this.

[ciaranevans]

Ah I see what you mean. I'll try give it a go tomorrow (finished work for today as I'm in the UK)

I should probably know this, but what's a easy way to make a dataset just based on bounds? I'm always just loading and transforming data rather than making it from scratch! 😓

[ciaranevans]

@sgillies been giving this a bit of a go.

Steps:

• Load in all tiles to 'merge'
• Find top, bottom, right and left (right goes +180 rather than wrapping to -180)
• Write a empty tiff out with those bounds that is the total size of the merged tiles (height + width), also has the same CRS as those tiles
• Load empty tiff, load tiles
• Reproject each tile into the big empty one
  Code for reprojecting below, where dst is the empty tiff and hgt(1/2) are the tiles.

dst_arr = dst.read()
kwargs = {'CENTER_LONG': 180}
warp.reproject(source=hgt1.read(),
                        destination=dst_arr,
                        src_transform=hgt1.transform,
                        src_crs=CRS.from_dict(init='epsg:4326'),
                        dst_transform=dst.transform,
                        dst_crs=CRS.from_dict(init='epsg:4326'),
                        resample=enums.Resampling.bilinear,
                        kwargs=kwargs)
show(dst_arr)

This is fine for the tile to the left of the AM:

But as soon as I move onto a tile the other side:

dst_arr = dst.read()
kwargs = {'CENTER_LONG': 180}
warp.reproject(source=hgt1.read(),
                        destination=dst_arr,
                        src_transform=hgt1.transform,
                        src_crs=CRS.from_dict(init='epsg:4326'),
                        dst_transform=dst.transform,
                        dst_crs=CRS.from_dict(init='epsg:4326'),
                        resample=enums.Resampling.bilinear,
                        kwargs=kwargs)
warp.reproject(source=hgt2.read(),
                        destination=dst_arr,
                        src_transform=hgt2.transform,
                        src_crs=CRS.from_dict(init='epsg:4326'),
                        dst_transform=dst.transform,
                        dst_crs=CRS.from_dict(init='epsg:4326'),
                        resample=enums.Resampling.bilinear,
                        kwargs=kwargs)
show(dst_arr)

I get an empty dataset:

I've tried with just doing the right hand side of the AM tiles and get the same thing

Any thoughts? 😄

[ciaranevans]

Not sure if this is the order of things you suggested to try!

[ciaranevans]

@sgillies I think I got it!

Annoyingly simple in the end!

Looks like a blob appeared on my AM!

Zoomed in without the basemap as my QGIS cries

Steps to reproduce

• Load SRTM HGT files
• For each tile:
  
  
  
  • Check if the left bound is negative
  
  
  • If negative, find the distance between 180 and the absolute value of the left bound, then set left to 180 + distance
  
  
  • Get top bound
  
  
  • Get old affine
  
  
  • Create a new affine with parameters (old.a, old.b, left, old.d, old.e, top_bound)
  
  
  • Write out as a tiff with new affine
  
  
• Load each tile into an array
• Use merge as normal

I'm gonna assume that this will happily reproject then into the MGRS CRS required, but it's late so I will 100% confirm that tomorrow.

[ciaranevans]

Sad times. For some reason QGIS reprojection works, see:
Merged tile in EPSG:32701

However, when I do the same with the following code in rasterio:

with rasterio.open('potentially_merged_2.tif', 'r') as merged:

    merged_full = merged.read()
    merged_full = np.clip(merged_full, 0, 400) # So colour ramp isn't awful
    show(merged_full)

    aff, w, h = warp.calculate_default_transform(
        src_crs=CRS.from_dict(init='epsg:4326'),
        dst_crs=CRS.from_dict(init=TILECRS),
        width=merged.meta['width'],
        height=merged.meta['height'],
        left=merged.bounds.left,
        right=merged.bounds.right,
        top=merged.bounds.top,
        bottom=merged.bounds.bottom,
        resolution=10
    )

    dest_arr = np.zeros((1, h, w), dtype='int16')

    warp.reproject(
        source=merged.read(),
        destination=dest_arr,
        src_transform=merged.transform,
        src_crs=CRS.from_dict(init='epsg:4326'),
        dst_transform=aff,
        dst_crs=CRS.from_dict(init=TILECRS),
        resampling=enums.Resampling.bilinear,
        kwargs={'CENTER_LONG':180}
    )
    show(dest_arr)

First tile is merged, second is reprojected

It's almost as if the width returned by calculate_default_transform() is too small, if I swap:

dest_arr = np.zeros((1, h, w), dtype='int16')

With:

dest_arr = np.zeros((1, h, h), dtype='int16')

And use h wherever I referred to w, I get:

My assumption is that calculate_default_transform() clips at the edge of the CRS...

@sgillies is their an alternative to figuring out the affine, height and width?

[ciaranevans]

Ah. Read some more into the docs for calculate_default_transform() and only noticed the Notes section. Also noticed that CHECK_WITH_INVERT_PROJ is defaulted to True.

Tried with:
with rasterio.Env(CHECK_WITH_INVERT_PROJ=False): wrapping it all and:

Which doesn't mess with the height or width.

Good to know that exists, this whole Anti Meridian issue is definitely a learning experience 🤣

[ciaranevans]

Summary

So, I have finally, finally got a perfect tile. We found that calculate_default_transform() didn't like the AM as our tiles were being shifted hundreds of meters off from the ground truth. I think this was due to the fact that these tiles stretch over two UTM Zones (Sad trombone) 🎺

We've got a work around though:

1. Merge tiles as per https://github.com/mapbox/rasterio/issues/1861#issuecomment-582540168
2. Get WGS84 coords for the BBox of the Sentinel 2 tile (Making sure to push Westings over 180 rather than following normal -180 to 0 convention)
3. Mask (with crop=True) the merged tiles with the BBox
4. Get the Affine for the Sentinel tile (We did this by loading one of the 10m .jp2 tiles from https://registry.opendata.aws/sentinel-2/)
5. Create a numpy array of zeros, dimensions=(1, s2_tile_height, s2_tile_width)
6. Reproject the cropped WGS84 image into the Sentinel tiles CRS (with kwargs={'CENTER_LONG:180})
7. Save output

Outputs

Sentinel 2 Scene with Pink border indicating MGRS bounds AND the Anti Meridian

Merged SRTM Tiles - Cropped to fit the shape of the MGRS bounds (With Hillshade)

Comparison of SRTM and Sentinel True Colour Image:

Sorry for the spamming of this issue @sgillies! Hopefully anyone in a similar situation will find this all handy. Not sure what the artefacting is at the bottom of our SRTM tile, but in our downstream pipelines that will get replaced by adjacent images.

❤️ You Rasterio

[geotheory]

Most interesting, but just WRT the title the International Date Line is not the same thing as 180 meridian.

[ciaranevans]

Most interesting, but just WRT the title the International Date Line is not the same thing as 180 meridian.

Yes good shout, sorry we'd been using it colloquially at work but have changed now!

[sgillies]

Thanks for summarizing @ciaranevans. I'm going to remove the bug label as merge can't merge at the edges of a reference system and the solution involves reprojecting the data.

[ciaranevans]

No worries @sgillies thanks for the input through it!

[geotheory]

Yes good shout, sorry we'd been using it colloquially at work but have changed now!

Ah good. I can stop being offended now.