Cartopy CRS definition for NZTM2000

cartopy-nztm.py

import cartopy
import shapely.geometry as sgeom

class NZTM(cartopy.crs.TransverseMercator):
    def __init__(self):
        '''
        https://epsg.io/2193
        https://www.linz.govt.nz/data/geodetic-system/datums-projections-and-heights/geodetic-datums/new-zealand-geodetic-datum-2000-nzgd2000
        '''
        inverse_flattening = 298.257222101
        semimajor_axis = 6378137
        semiminor_axis = -((inverse_flattening*semimajor_axis)-semimajor_axis)
        globe = cartopy._crs.Globe(
            semimajor_axis=semimajor_axis,
            semiminor_axis=semiminor_axis
        )
        super(NZTM, self).__init__(central_longitude=173,
                                   central_latitude=0,
                                   scale_factor=0.9996,
                                   false_easting=1600000,
                                   false_northing=10000000,
                                   globe=globe)

    @property
    def boundary(self):
        x0, x1 = self.x_limits
        y0, y1 = self.y_limits
        return sgeom.LineString([(x0, y0), (x0, y1),
                                 (x1, y1), (x1, y0),
                                 (x0, y0)])

    @property
    def x_limits(self):
        return (827933.23, 3195373.59)

    @property
    def y_limits(self):
        return (3729820.29, 7039943.58)

For example, in a Flask application, to return a PNG of NZ in NZTM to the client at http://host:port/map:

@app.route('/map')
def getmap():
    ax = matplotlib.pyplot.axes(projection=NZTM())
    ax.add_feature(cartopy.feature.LAND)
    ax.coastlines('10m', alpha=0.2)
    x0, y0 = 165.893555, -47.696163
    x1, y1 = 180.0, -34.016242
    ax.set_extent([x0, x1, y0, y1])
    plt.title('New Zealand Transverse Mercator (2000) with Cartopy')
    output = io.BytesIO()
    matplotlib.backends.backend_agg.FigureCanvasAgg(ax.figure).print_png(output)
    return Response(output.getvalue(), mimetype='image/png')