Mapping proj4 to cartopy CRS

[pelson]

Consider the feasibility of mapping proj4 definitions to cartopy classes / keywords.

This may be hard, but would yield several benefits, including the ability to get users new to cartopy (but not proj4) up to speed with the class hierarchy.

Another benefit is the following file, along with other standard osgeo projects which talk in proj4 definitions: https://github.com/matplotlib/basemap/blob/master/lib/mpl_toolkits/basemap/data/epsg

[rhattersley]

Seems feasible to map proj.4 definitions to Cartopy instances, but it won't always give the expected result[1] as proj.4 definitions don't include all the information that Cartopy uses (e.g. domain of use). WKT should be a richer source.

[1] - The proj.4 definition of EPSG:27700 (British National Grid) doesn't include the domain of validity, so a generic conversion would just result in a transverse Mercator.

[rhattersley]

In fact I've been wondering if Cartopy's class structure would be well served by mimicking that of WKT.

[pelson]

#275 adds the ability to specify a datum/ellipse.

[snorfalorpagus]

I came across a need for this when trying to display an image parsed from rasterio using cartopy. I managed to shoehorn the relevant information from the rasterio image metadata into the correct cartopy CRS object. It seems a lot of images could be handled this way with the use of a helper function of some sort.

What is it in particular that makes this feature hard to implement?

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
import cartopy.crs as ccrs
import rasterio
import rasterio.features
import matplotlib.pyplot as plt

if __name__ == '__main__':
    filename = 'RGB.byte.tif' # from the rasterio test data
    with rasterio.open(filename, 'r') as src:
        meta = src.meta
        width = meta['width']
        height = meta['height']
        count = meta['count']
        dtype = meta['dtype']

        # allocate memory for image
        im = np.empty([height, width, count], dtype)

        # read image into memory
        for band in range(0, meta['count']):
            im[:,:,band] = src.read_band(band+1)

        # calculate extent of raster
        t = meta['transform']
        xmin = t[0]
        xmax = t[0] + t[1]*meta['width']
        ymin = t[3] + t[5]*meta['height']
        ymax = t[3]

        # define cartopy crs for the raster, based in rasterio metadata
        crs = ccrs.UTM(meta['crs']['zone'])

        # create figure
        ax = plt.axes(projection=crs)
        plt.title('RGB.byte.tif')
        ax.set_xmargin(0.05)
        ax.set_ymargin(0.10)

        # plot raster
        plt.imshow(im, origin='upper', extent=[xmin, xmax, ymin, ymax], transform=crs)

        # plot coastlines
        ax.coastlines(resolution='10m', color='red', linewidth=1)

        plt.show()

[pelson]

What is it in particular that makes this feature hard to implement?

You've got a projection which is easy to map to a cartopy crs (crs = ccrs.UTM(meta['crs']['zone'])).
To solve this generally, there are a lot of parameters which make up a complex Proj4 definition, and they need to be mapped to the appropriate keyword on the appropriate class (additionally some of the parameters need to be manipulated into a more sane form). All in all solving the problem generally is hard - but certainly doable.

In fact, I have a branch which proves the concept - it just needs a bit of work to polish it off and it would make a nice addition to cartopy (I'm taking a long flight in a fortnight and planning on working on this and another major feature for matplotlib with a bit of luck).

I came across a need for this when trying to display an image parsed from rasterio using cartopy.

One of the biggest benefits of implementing the proj4 <-> cartopy CRS interface is that we can instantly hook in to some awesome projects (such as rasterio) and get a whole host of very powerful functionality for free.

Great example by the way - something I'd love to see in Cartopy's gallery, maybe before even having the proj4 mappings... is that something you'd be interested in submitting?

[snorfalorpagus]

I'd be happy to help implement this, even if it's just doing some testing.

Re: an example for the gallery, I think it could fit well - maybe add some simple raster processing to illustrate why you might want to do this?

[pelson]

maybe add some simple raster processing to illustrate why you might want to do this?

Perhaps - but to be honest a simple example is a great start. I've got my eye on publishing a series of worked examples (in the form of IPython notebooks) which would be a much better fit for an indepth analyses (in my head some interesting image processing would be perfect there).

As I say, simple is good to start with, but either way, your contribution would be most welcome.

Cheers!

[marqh]

I think this is a good time to be thinking about interfacing to well known text.

The OGC standards working group on well known text recently published their revamped standard, which is being adopted by OGC and is being proposed as an ISO standard (i'll provide a public URI when I find one).
This should provide a stable platform to code against.

using GDAL (if it can still parse WKT) to produce proj4 and extending cartopy so it can manage a
proj4 -> ccrs
process sounds like a neat way to interface cartopy to proj4 and WKT.

rhattersley's

In fact I've been wondering if Cartopy's class structure would be well served by mimicking that of WKT.

sounds like an interesting approach, with this in mind

[naught101]

Is there currently a recommended way to get from a GDAL spatial reference to a cartopy CRS? is WKT the best way?