TL;DR

Python example here.

NZTM to Web Mercator

There are a number of times when I am building a web map when I need to warp a GTiff across the anti-meridian. New Zealand sits at the edge of this line, with a portion of it on the other side, The Chatham Islands. Doing a direct gdalwarp from NZTM (2193) to Web Mercator (3857) will “split” the file with the portion across the anti-meridian being placed to the other side of the projection.

Doing a straight warp from NZTM to Web Mercator:

gdalwarp -s_srs EPSG:2193 -t_srs:3857 in_tif.tif webmer_out_tif.tif

produces the following:

A number of people have devised methods to get around this, but I wanted to post one way I found that could be useful, albeit verbose. This is in by no means the best method; however, I have found it useful as a solution when I have a request to build a visualization using Web Mercator.

The short description is this:

The goal is to get the extents in Web Mercator for your GTiff that cross the anti-meridian. These extents can be passed to gdal.Warp as the output extents in your warp command

A Github repo with an example script can be found here. The code is intentionally verbose to help keep track of what is happening along the way.

The example uses Python and a number of available libraries, including:

• geopandas
• gdal
• rasterio
• fiona

The real star of the method is fiona. Fiona has a method in the command that allows reprojection across the anti-meridian: antimeridian_cutting=True. When using this, fiona will output a multi-polygon in two pieces similar to the image above. Using the individual extents from each polygon in the multi-polygon, the extents can be leveraged to reconstruct a polygon across the anti-meridian. The reconstructed polygon extents can then be used as the output extents of you gdalwarp operation.

A little more in depth:

• Find extent of input file
• Reproject this extent to WGS 84 using fiona with antimeridian_cutting=True. This creates a multi-polygon
• Find extent of each multi-poly
• Shift left polygon over anti-meridian, next to right polygon
• Build extent of polygons combined
• Reproject WGS polygon to Web Mercator
• Shift maxx of Web Mercator polygon across the anti-meridian
• Find extent of new Webmercator polygon
• Use new extent as output bounds in gdal.Warp process

Again, here is link to the repo with an example.

The example in the repo uses data from National Institute of Water and Atmospheric Research (NIWA) in New Zealand. You can get a download of the example data via Koordinates here.

Normally, I would try to build the code in the blog post, but the example too large. Also, I usually try to break down all the import highlights in the code; however, there is just too much to cover in the method. Hopefully, this will get you a little further along.