New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Mapping proj4 to cartopy CRS #153

Open
pelson opened this Issue Nov 15, 2012 · 8 comments

Comments

Projects
None yet
4 participants
@pelson
Member

pelson commented Nov 15, 2012

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

This comment has been minimized.

Show comment
Hide comment
@rhattersley

rhattersley Nov 15, 2012

Member

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.

Member

rhattersley commented Nov 15, 2012

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

This comment has been minimized.

Show comment
Hide comment
@rhattersley

rhattersley Nov 15, 2012

Member

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

Member

rhattersley commented Nov 15, 2012

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

@pelson pelson referenced this issue Feb 7, 2013

Closed

wms #114

@pelson

This comment has been minimized.

Show comment
Hide comment
@pelson

pelson Jun 3, 2013

Member

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

Member

pelson commented Jun 3, 2013

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

@snorfalorpagus

This comment has been minimized.

Show comment
Hide comment
@snorfalorpagus

snorfalorpagus Jun 24, 2014

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()

figure_1

snorfalorpagus commented Jun 24, 2014

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()

figure_1

@pelson

This comment has been minimized.

Show comment
Hide comment
@pelson

pelson Jun 25, 2014

Member

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?

Member

pelson commented Jun 25, 2014

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

This comment has been minimized.

Show comment
Hide comment
@snorfalorpagus

snorfalorpagus Jun 25, 2014

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?

snorfalorpagus commented Jun 25, 2014

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

This comment has been minimized.

Show comment
Hide comment
@pelson

pelson Jun 26, 2014

Member

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!

Member

pelson commented Jun 26, 2014

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

This comment has been minimized.

Show comment
Hide comment
@marqh

marqh Sep 30, 2014

Member

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

Member

marqh commented Sep 30, 2014

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment