Add support for cartopy projections with geographic data #1966

munkm · 2018-08-10T23:03:41Z

PR Summary

This PR adds in support for projecting geographic data sources on various map types. As it stands, the code in this PR assumes that geographic data is read in as PlateCarree type and can be projected onto the maps listed here: https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.

PR Checklist

Code passes flake8 checker
New features are documented, with docstrings and narrative docs
Adds a test for any bugs fixed. Adds tests for new features.

matthewturk · 2018-08-16T20:10:50Z

Cartopy projections for tomography stuff

munkm · 2018-08-16T20:20:14Z

Update: I'm still finishing up some stuff for this PR. It probably shouldn't be merged until after the cf-radial branch has been PR'd and merged since I've been periodically merging the work there into this branch.

munkm · 2018-08-22T22:58:09Z

Ok, I think this is pretty close to a working state.

This branch:

Changes the axes for a geographic coordinate system-based dataset to be an instance of cartopy.mpl.geoaxes.GeoAxes. These axes can then be annotated as you would with cartopy.
The default transform and projection type is assumed to be PlateCarree
includes a new function set_mpl_projection() which can be used to set the projection on a dataset that has a detected geographic coordinate system

a default slice plot with this is PlateCarree

ds = yt.load_uniform_grid(data, sizes, 1.0, geometry=("geographic", dims),
                          bbox=bbox)
p = yt.SlicePlot(ds, "altitude", 'AIRDENS')
p.show()

This can be annotated like any other plot with GeoAxes to give contextual information.

p._setup_plots()
p.plots['AIRDENS'].axes.set_global()
p.plots['AIRDENS'].axes.coastlines()
p.show()

now set_mpl_projection() can be used to transform and project this data onto new axes

a few examples:

p.set_mpl_projection('AlbersEqualArea')
p._setup_plots()
p.plots['AIRDENS'].axes.set_global()
p.plots['AIRDENS'].axes.coastlines()
p.show()

and

p.set_mpl_projection('InterruptedGoodeHomolosine')
p._setup_plots()
p.plots['AIRDENS'].axes.set_global()
p.plots['AIRDENS'].axes.coastlines()
p.show()

ngoldbaum · 2018-08-22T23:01:16Z

ping @dopplershift and @pelson you guys might be interested in this :)

ngoldbaum · 2018-08-22T23:02:50Z

@zssherman how do you want to deal with upstreaming the cf_radial frontend? That should probably happen separately from this.

ngoldbaum · 2018-08-22T23:04:31Z

I think we probably want to add a way to call set_global() and coastlines() on the plot object itself, but that can come separately, I guess via the plot callback that @zssherman was working on.

This is really, really awesome.

zssherman · 2018-08-23T00:00:06Z

@ngoldbaum I can do it separately if you would like? Sadly I've been out sick, so it might take a day or two. I do have coastlines etc in the plot callback. And @munkm Great work!

munkm · 2018-08-23T02:24:28Z

Thanks @zssherman and @ngoldbaum! It means a lot coming from the both of you!

munkm · 2018-08-23T02:32:18Z

So I meant to bring up a few issues I have with this that I'd like to ask about before I take off the 'WIP' label, one of which @ngoldbaum already brought up!

I've been periodically merging in @zssherman 's cf_radial branch. Should I wait to take off the WIP label until that's merged?
This work is all technically independent from the cf_radialbranch and could probably be merged separately, but I didn't want to do that unless we agree that it's the way to proceed.
I have the unit tests comparing the plot with GenericImageTest and I also do some tests with typechecking that the projection is set in our plots. The unit test usingGenericImageTest might be redundant but I though it might be useful to compare plots as well as the projection types. Let me know if these tests should be restructured.
I have added documentation, flake8 checks, and unit tests for the features I've added, so I checked all of those off in the PR checklist, but the CI won't pass because I haven't totally updated everything from the cf_radial branch.
I'm also happy to work with @zssherman on the plot callback annotations, but I'd like some advice on how to proceed with this PR first (i.e. how should things get merged? Should the callbacks be in this PR? Should the cf_radial get merged as is and I add the features for the callback in this PR for projections? Which branch should have the annotations with the transform features?)

munkm · 2018-08-23T03:08:12Z

Also I'm sorry to hear that you're out sick @zssherman! I hope you feel better soon. 🙂

munkm · 2018-08-24T19:41:42Z

Ok I rebased so that the cf-radial branch can be submitted as a separate PR by @zssherman. I think this should separate out the projection work from the cf-radial frontend.

matthewturk

Looks good to me. Only minor comments.

It is worth noting that down the road we will likely want to shift the spherical coordinate handler to this, too, even though it won't necessarily need coastlines.

What do you think about making a note if the data is geographic and the projection is not set, and saying, hey, you might want to set the projection somehow.

matthewturk · 2018-08-24T21:34:51Z

doc/source/visualizing/geographic_projections_and_transforms.rst

+onto a representation of that sphere flattened into 2d space. There exist a
+number of projection types, which can be found in the `the cartopy
+documentation <https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html>`_.
+yt now supports these projection types for geographically loaded data.


Can you change this to be:

`yt` now supports using these projection types, through cartopy, for geopgraphically loaded data.

or something that indicates more directly we're using cartopy?

Ah, great point! Cartopy is definitely doing all of the transformation behind-the-scenes so it should be highlighted.

matthewturk · 2018-08-24T21:35:15Z

doc/source/visualizing/geographic_projections_and_transforms.rst

+number of projection types, which can be found in the `the cartopy
+documentation <https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html>`_.
+yt now supports these projection types for geographically loaded data.
+Underlying data is assumed to have an underlying projection type of `PlateCarree`. If


Make this a link to the docs for PlateCaree, or at least note that PlateCaree means that it's lat/lon grid?

matthewturk · 2018-08-24T21:35:31Z

doc/source/visualizing/geographic_projections_and_transforms.rst

+request on the yt github page for this feature.
+
+It should be noted that
+these projections are not the same as yt's ProjectionPlot. For more information


matthewturk · 2018-08-24T21:36:33Z

doc/source/visualizing/geographic_projections_and_transforms.rst

+
+.. code-block:: python
+
+    p.set_mpl_projection('Robinson')


Is it possible to do this in the instantiation of SlicePlot itself, with an argument?

Right now I don't think so, but I could add that if it'd be useful.

I think it would be good, so that we can enable it to be done as the first result, rather than as the second iteration.

Ok, I've added in that functionality and updated the test and the documentation to reflect this additional functionality. We should now be able to pass in geo_projection= into SlicePlot and the initial plot should have a projection other than PlateCarree.

munkm · 2018-11-12T19:28:26Z

@ngoldbaum and I just talked about updating the docs on this PR to ensure that users wanting to use this feature have cartopy installed to avoid potential issues with the GEOS library in shapely / cartopy. We've decided that the following improvements would be helpful:

adding in additional clarification on installing cartopy with various package management options in the geographic_projections_and_transforms.rst file. This subsection on installation will be referenced in the yt general installation instructions and also in the geographic transforms cookbook.
adding additional clarification to the warnings that'll be generated if somebody doesn't have cartopy installed. If the user is using conda this message will be relatively short. If not the error will be a bit longer and detail needing to install GEOS and proj4 with their system package manager and build shapely from source.

Nathan, I think that covers everything. Did I miss anything?

ngoldbaum · 2018-11-12T19:30:37Z

Yup and just because it was annoying to figure out, the pip install command that should work on macs is:

pip install --no-binary :all: shapely cartopy

matthewturk · 2018-11-12T21:28:33Z

Sounds good. Maybe we could out an error message if a projection is chosen that is not supported without cartopy.

matthewturk · 2018-11-13T20:17:52Z

I'm running into some issues when plotting data that doesn't cover the whole domain.

I'm using this dataset: http://ds.iris.edu/ds/products/emc-dna13/ the binary version of which is at http://ds.iris.edu/files/products/emc/data/DNA13/DNA13_percent.nc . ( Porritt, R. W., R.M. Allen, and F. F. Pollitz. 2014. “Seismic imaging east of the Rocky Mountains with USArray.” Earth Planet. Sci. Lett. https://doi.org/10.1016/j.epsl.2013.10.034. )

Here is a gist that loads into yt:

https://gist.github.com/MatthewTurk/57360f40433817ecee4ef76bb817d970

It hangs inside cartopy's img_transform.py in the call to kdtree.query. I inserted some print-debuggers that suggest the target_xyz variable has a maximum value of inf along the axis=0 in the first and second locations.

I'm not entirely certain why this is happening. Cartopy seems to be expecting some distances that are inf, but not before the kdtree query call. I'm struggling to make this an outside-of-yt reproducible example, which is why I'm posting it here. If I am able to fix it on our end I will post an update in the code.

matthewturk · 2018-11-13T20:39:10Z

OK, I was able to reduce it to something that doesn't use yt; now I'm trying to determine how we might avoid it in the future.

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

data = np.random.random((256, 256))

ax = plt.axes(projection=ccrs.Mollweide())
extent = [-27.0, 27.0, -13.5, 13.5]
ax.set_extent(extent)
ax.imshow(data, extent = extent, transform = ccrs.PlateCarree())

plt.savefig('example.png')

matthewturk · 2018-11-13T21:02:29Z

I submitted an issue upstream here: SciTools/cartopy#1184

pelson · 2018-11-14T05:22:44Z

Hope my clarification in the issue above was helpful. In summary, it is just taking a really long time for cartopy to re-project the data. v0.17 (release due this week) of cartopy will support pykdtree, and when installed will make the transformation an order of magnitude quicker. Hopefully that solves the issue in yt for you @munkm.

matthewturk · 2018-11-14T14:21:41Z

Seems that my "hotfixes" for the extent-setting have caused secondary problems; I'm going to revert them and then we can figure out how to move forward with either the forthcoming new release of cartopy that @pelson noted or if I can find an alternate implementation.

brittonsmith · 2018-11-27T18:07:55Z

doc/source/visualizing/geographic_projections_and_transforms.rst

+^^^^^^^^^^^^^^^^^^^^^^^
+
+As mentioned above, the default data transform is assumed to be of `PlateCarree
+<https://scitools.org.uk/cartopy/docs/latest/crs/projections.html#platecarree>`_,


This is a very minor note. There are multiple places in this doc file with a link where the link text is "PlateCarree." This will cause sphinx to think a label has multiple definition and issue a warning. This can be fixed by making the trailing underscore into two underscores, e.g.:

`PlateCarree <https://scitools.org.uk/cartopy/docs/latest/crs/projections.html#platecarree>`__

Ah, thanks for that @brittonsmith! I really appreciate the explanation too. It should be fixed now.

brittonsmith

This looks great to me!

brittonsmith · 2018-11-27T18:08:39Z

doc/source/visualizing/geographic_projections_and_transforms.rst

+``yt`` now supports these projection 
+types for geographically loaded data.
+Underlying data is assumed to have a transform of `PlateCarree
+<https://scitools.org.uk/cartopy/docs/latest/crs/projections.html#platecarree>`_,


There's the other one.

ngoldbaum · 2018-11-27T20:19:47Z

@munkm along with britton's suggestion, maybe we should try testing with cartopy 0.17.0 now that that has been released. Supposedly that will fix the perf issues that Matt ran into. If updating cartopy causes issues we can back that out and go in and fix those issues in a separate PR.

ngoldbaum · 2018-11-27T20:19:59Z

BTW with Britton approving I feel comfortable merging this.

matthewturk and others added 15 commits August 24, 2018 14:27

Experimentally adding mpl_proj argument

5b14e6a

Changing to default mpl_proj

2b9302a

WIP: geo projection plotting

4e0e8c4

getting geo projections to pass to axes

2bfecbc

add in data transform functionality to geo data

a963d99

modify geographic pixelization from aitoff to cartesian

83be059

move cartopy helper functions to standalone file

a1e5d60

add in remaining cartopy projections

c205a41

add docstrings to new funcs

e4dfbed

remove unused import

ccb5872

fix misspelled transform names

9fb4c69

add unittests for new feature

4fd5417

add docs for geographic projections in visualization section

935d742

add cookbook entry and notebook for geo projections

8395d8b

fix flake8 stuff

ebe7eb0

munkm force-pushed the yt_mpl_proj branch from 7a776b9 to ebe7eb0 Compare August 24, 2018 19:39

munkm changed the title ~~[WIP] Add support for cartopy projections with geographic data~~ Add support for cartopy projections with geographic data Aug 24, 2018

matthewturk approved these changes Aug 24, 2018

View reviewed changes

add in cartopy docs clarification

b2c84fd

Nathan Goldbaum added 5 commits November 9, 2018 14:37

don't use the preinstall numpy from homebrew

c20ad8c

uninstall things that need numpy too

5b4e648

avoid permissions error

0f430d0

one more try to fix permissions issue

e921e23

skip blacklisted cartopy projections

0e4c59d

skip UTM in answer tests as well

7e0ae14

ngoldbaum pushed a commit to ngoldbaum/yt that referenced this pull request Nov 13, 2018

more fixes from yt-project#1966

9235ecf

ngoldbaum mentioned this pull request Nov 13, 2018

Travis build improvements #2093

Merged

matthewturk added 2 commits November 13, 2018 14:46

Setting extent before imshow

b8f15b0

Only call set_extent if it exists

e1d6596

matthewturk and others added 4 commits November 16, 2018 10:00

semi-revert set_extent

c0341eb

Merge branch 'master' into yt_mpl_proj

f294a9e

add cartopy installation directions to projection docs

e30770f

return system-dependent error message if cartopy is not installed

5b0f68d

brittonsmith reviewed Nov 27, 2018

View reviewed changes

brittonsmith approved these changes Nov 27, 2018

View reviewed changes

munkm added 3 commits November 27, 2018 15:03

change trailing underscores on PlateCarree links in docs

d130d65

update cartopy version to 0.17.0 for tests

f6f1a19

update geo answer test

6e8ba0a

ngoldbaum merged commit ba2b2ef into yt-project:master Nov 28, 2018

Add support for cartopy projections with geographic data #1966

Add support for cartopy projections with geographic data #1966

Uh oh!

Conversation

munkm commented Aug 10, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

PR Summary

PR Checklist

Uh oh!

matthewturk commented Aug 16, 2018

Uh oh!

munkm commented Aug 16, 2018

Uh oh!

munkm commented Aug 22, 2018

Uh oh!

ngoldbaum commented Aug 22, 2018

Uh oh!

ngoldbaum commented Aug 22, 2018

Uh oh!

ngoldbaum commented Aug 22, 2018

Uh oh!

zssherman commented Aug 23, 2018 via email • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

munkm commented Aug 23, 2018

Uh oh!

munkm commented Aug 23, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

munkm commented Aug 23, 2018

Uh oh!

munkm commented Aug 24, 2018

Uh oh!

matthewturk left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

munkm commented Nov 12, 2018

Uh oh!

ngoldbaum commented Nov 12, 2018

Uh oh!

matthewturk commented Nov 12, 2018

Uh oh!

matthewturk commented Nov 13, 2018

Uh oh!

matthewturk commented Nov 13, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

matthewturk commented Nov 13, 2018

Uh oh!

pelson commented Nov 14, 2018

Uh oh!

matthewturk commented Nov 14, 2018

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

brittonsmith left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ngoldbaum commented Nov 27, 2018

Uh oh!

munkm commented Aug 10, 2018 •

edited

Loading

zssherman commented Aug 23, 2018 via email •

edited

Loading

munkm commented Aug 23, 2018 •

edited

Loading

matthewturk commented Nov 13, 2018 •

edited

Loading