Skip to content
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

Fix bugs and use real data in the plotting docs #251

Merged
Merged
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions docs/source/geo_def.rst
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@ classes must be in degrees. Additionally, longitudes must be in the
AreaDefinition
--------------

.. _area-definitions:

An :class:`~pyresample.geometry.AreaDefinition`, or ``area``, is the primary
way of specifying a uniformly spaced geographic region in pyresample. It is
also one of the only geometry objects that understands geographic projections.
Expand Down
165 changes: 96 additions & 69 deletions docs/source/plot.rst
Original file line number Diff line number Diff line change
Expand Up @@ -22,30 +22,48 @@ supported by Satpy_. All you need are a set of geo-referenced values

First we read in the data with Satpy_:

>>> from satpy.scene import Scene
>>> from glob import glob
>>> SCENE_FILES = glob("./GW1AM2_20191122????_156*h5")
>>> scn = Scene(reader='amsr2_l1b', filenames=SCENE_FILES)
>>> scn.load(["btemp_36.5v"])
>>> lons, lats = scn["btemp_36.5v"].area.get_lonlats()
>>> tb37v = scn["btemp_36.5v"].data.compute()
>>> from satpy.scene import Scene # doctest: +SKIP
>>> from glob import glob # doctest: +SKIP
>>> SCENE_FILES = glob("./GW1AM2_20191122????_156*h5") # doctest: +SKIP
>>> scn = Scene(reader='amsr2_l1b', filenames=SCENE_FILES) # doctest: +SKIP
>>> scn.load(["btemp_36.5v"]) # doctest: +SKIP
>>> lons, lats = scn["btemp_36.5v"].area.get_lonlats() # doctest: +SKIP
>>> tb37v = scn["btemp_36.5v"].data.compute() # doctest: +SKIP

Data for this example can be downloaded from zenodo_.

If you have your own data, or just want to see that the example code here runs, you can
set the three arrays :code:`lons`, :code:`lats` and :code:`tb37v` accordingly, e.g.:

.. doctest::
:hide:
>>> from pyresample.geometry import AreaDefinition
>>> area_id = 'ease_sh'
>>> description = 'Antarctic EASE grid'
>>> proj_id = 'ease_sh'
>>> projection = {'proj': 'laea', 'lat_0': -90, 'lon_0': 0, 'a': 6371228.0, 'units': 'm'}
>>> width = 425
>>> height = 425
>>> area_extent = (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)
>>> area_def = AreaDefinition(area_id, description, proj_id, projection,
... width, height, area_extent)

>>> import numpy as np
>>> from pyresample import load_area, save_quicklook, SwathDefinition
>>> from pyresample.kd_tree import resample_nearest
>>> lons = np.zeros(1000)
>>> lats = np.arange(-80, -90, -0.01)
>>> tb37v = np.arange(1000)
>>> area_def = load_area('areas.yaml', 'ease_sh')

But here we go on with the loaded AMSR-2 data. Make sure you have an :code:`areas.yaml`
file that defines the :code:`ease_sh` area. Or see
:ref:`the area definition section<area-definitions>` on how to define one.

>>> from pyresample import load_area, save_quicklook, SwathDefinition
>>> from pyresample.kd_tree import resample_nearest
>>> area_def = load_area('areas.yaml', 'ease_sh') # doctest: +SKIP
>>> swath_def = SwathDefinition(lons, lats)
>>> result = resample_nearest(swath_def, tb37v, area_def,
... radius_of_influence=20000, fill_value=None)
>>> save_quicklook('tb37v_quick.png', area_def, result, label='Tb 37v (K)')

Assuming **lons**, **lats** and **tb37v** are initialized with real data (as in
the above Satpy_ example) the result might look something like this:

Expand Down Expand Up @@ -78,18 +96,23 @@ Assuming the file **areas.yaml** has the following area definition:
**Example usage:**

.. doctest::
:hide:
>>> from pyresample.geometry import AreaDefinition
>>> area_id = 'pc_world'
>>> description = 'Plate Carree world map'
>>> proj_id = 'eqc'
>>> projection = {'proj': 'eqc', 'lat_0': -40, 'lon_0': 40, 'a': 6370997.0, 'units': 'm'}
>>> width = 640
>>> height = 480
>>> area_extent = (-20037508.34, -10018754.17, 20037508.34, 10018754.17)
>>> area_def = AreaDefinition(area_id, description, proj_id, projection,
... width, height, area_extent)

>>> import numpy as np
>>> from pyresample import load_area, save_quicklook, SwathDefinition
>>> from pyresample.kd_tree import resample_nearest
>>> lons = np.zeros(1000)
>>> lats = np.arange(-80, -90, -0.01)
>>> tb37v = np.arange(1000)
>>> area_def = load_area('areas.yaml', 'pc_world')
>>> swath_def = SwathDefinition(lons, lats)
>>> area_def = load_area('areas.yaml', 'pc_world') # doctest: +SKIP
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Isn't this area load just overwriting the area_def created in the lines above?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, right. But the stuff above is hidden, and does not show up on RTD. Therefore the additional line for RTD. The area is not defined and on Travis we don't have access to areas.yaml and certainly not in the home dir.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, pyresample doesn't have a builtin areas 🤦‍♂️ Got it.

>>> result = resample_nearest(swath_def, tb37v, area_def, radius_of_influence=20000, fill_value=None)
>>> save_quicklook('tb37v_pc.png', area_def, result, num_meridians=None, num_parallels=None, label='Tb 37v (K)')


Assuming **lons**, **lats** and **tb37v** are initialized with real data (like
above we use AMSR-2 data in this example) the result might look something like
this:
Expand Down Expand Up @@ -124,15 +147,22 @@ the following area definition for an ortho projection area:
**Example usage:**

.. doctest::

>>> import numpy as np
:hide:
>>> from pyresample import load_area, save_quicklook, SwathDefinition
>>> from pyresample.kd_tree import resample_nearest
>>> lons = np.zeros(1000)
>>> lats = np.arange(-80, -90, -0.01)
>>> tb37v = np.arange(1000)
>>> area_def = load_area('areas.yaml', 'ortho')
>>> from pyresample.geometry import AreaDefinition
>>> area_id = 'ortho'
>>> description = 'Ortho globe'
>>> proj_id = 'ortho'
>>> projection = {'proj': 'ortho', 'lat_0': -40, 'lon_0': 40, 'a': 6370997.0, 'units': 'm'}
>>> width = 640
>>> height = 480
>>> area_extent = (-10000000, -10000000, 10000000, 10000000)
>>> area_def = AreaDefinition(area_id, description, proj_id, projection,
... width, height, area_extent)

>>> swath_def = SwathDefinition(lons, lats)
>>> area_def = load_area('areas.yaml', 'ortho') # doctest: +SKIP
>>> result = resample_nearest(swath_def, tb37v, area_def, radius_of_influence=20000, fill_value=None)
>>> save_quicklook('tb37v_ortho.png', area_def, result, num_meridians=None, num_parallels=None, label='Tb 37v (K)')

Expand All @@ -153,33 +183,27 @@ converted to a Cartopy CRS as long as Cartopy can represent the
projection. Once an `AreaDefinition` is converted to a CRS object it can be
used like any other Cartopy CRS object.

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from pyresample import load_area, SwathDefinition
>>> from pyresample.kd_tree import resample_nearest
>>> from pyresample.geometry import AreaDefinition
>>> lons = np.zeros(1000)
>>> lats = np.arange(-80, -90, -0.01)
>>> tb37v = np.arange(1000)
>>> swath_def = SwathDefinition(lons, lats)
>>> import matplotlib.pyplot as plt # doctest: +SKIP
>>> area_id = 'alaska'
>>> description = 'Alaska Lambert Equal Area grid'
>>> proj_id = 'alaska'
>>> projection = {'proj': 'stere', 'lat_0': 62., 'lon_0': -152.5, 'ellps': 'WGS84', 'units': 'm'}
>>> width = 2019
>>> height = 1463
>>> area_extent = (-757214.993104, -485904.321517, 757214.993104, 611533.818622)
>>> from pyresample.geometry import AreaDefinition
>>> area_def = AreaDefinition(area_id, description, proj_id, projection,
... width, height, area_extent)
>>> result = resample_nearest(swath_def, tb37v, area_def,
... radius_of_influence=20000, fill_value=None)
>>> crs = area_def.to_cartopy_crs()
>>> ax = plt.axes(projection=crs)
>>> ax.coastlines()
>>> ax.set_global()
>>> plt.imshow(result, transform=crs, extent=crs.bounds, origin='upper')
>>> plt.colorbar()
>>> plt.savefig('amsr2_tb37v_cartopy.png')
>>> crs = area_def.to_cartopy_crs() # doctest: +SKIP
>>> ax = plt.axes(projection=crs) # doctest: +SKIP
>>> ax.coastlines() # doctest: +SKIP
>>> ax.set_global() # doctest: +SKIP
>>> plt.imshow(result, transform=crs, extent=crs.bounds, origin='upper') # doctest: +SKIP
>>> plt.colorbar() # doctest: +SKIP
>>> plt.savefig('amsr2_tb37v_cartopy.png') # doctest: +SKIP


Assuming **lons**, **lats**, and **i04_data** are initialized with real data
the result might look something like this:
Expand All @@ -199,21 +223,30 @@ AreaDefinition using the **plot.area_def2basemap(area_def, **kwargs)** function.

**Example usage:**

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from pyresample import load_area, save_quicklook, area_def2basemap, SwathDefinition
.. doctest::
:hide:
>>> from pyresample.kd_tree import resample_nearest
>>> lons = np.zeros(1000)
>>> lats = np.arange(-80, -90, -0.01)
>>> tb37v = np.arange(1000)
>>> area_def = load_area('areas.yaml', 'ease_sh')
>>> swath_def = SwathDefinition(lons, lats)
>>> from pyresample.geometry import AreaDefinition
>>> area_id = 'ease_sh'
>>> description = 'Antarctic EASE grid'
>>> proj_id = 'ease_sh'
>>> projection = {'proj': 'laea', 'lat_0': -90, 'lon_0': 0, 'a': 6371228.0, 'units': 'm'}
>>> width = 425
>>> height = 425
>>> area_extent = (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)
>>> area_def = AreaDefinition(area_id, description, proj_id, projection,
... width, height, area_extent)

>>> import matplotlib.pyplot as plt # doctest: +SKIP
>>> from pyresample import area_def2basemap # doctest: +SKIP
>>> area_def = load_area('areas.yaml', 'ease_sh') # doctest: +SKIP
>>> result = resample_nearest(swath_def, tb37v, area_def,
... radius_of_influence=20000, fill_value=None)
>>> bmap = area_def2basemap(area_def)
>>> bmng = bmap.bluemarble()
>>> col = bmap.imshow(result, origin='upper', cmap='RdBu_r')
>>> plt.savefig('tb37v_bmng.png', bbox_inches='tight')
>>> bmap = area_def2basemap(area_def) # doctest: +SKIP
>>> bmng = bmap.bluemarble() # doctest: +SKIP
>>> col = bmap.imshow(result, origin='upper', cmap='RdBu_r') # doctest: +SKIP
>>> plt.savefig('tb37v_bmng.png', bbox_inches='tight') # doctest: +SKIP


Assuming **lons**, **lats** and **tb37v** are initialized with real data as in
the previous examples the result might look something like this:
Expand All @@ -239,24 +272,18 @@ The above image can be generated using Cartopy instead by utilizing the method

**Example usage:**

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from pyresample import load_area, save_quicklook, area_def2basemap, SwathDefinition
>>> from pyresample.kd_tree import resample_nearest
>>> lons = np.zeros(1000)
>>> lats = np.arange(-80, -90, -0.01)
>>> tb37v = np.arange(1000)
>>> area_def = load_area('areas.yaml', 'ease_sh')
>>> swath_def = SwathDefinition(lons, lats)
>>> import matplotlib.pyplot as plt # doctest: +SKIP
>>> area_def = load_area('areas.yaml', 'ease_sh') # doctest: +SKIP
>>> result = resample_nearest(swath_def, tb37v, area_def,
... radius_of_influence=20000, fill_value=None)
>>> import matplotlib.pyplot as plt
>>> crs = area_def.to_cartopy_crs()
>>> ax = plt.axes(projection=crs)
>>> ax.background_img(name='BM')
>>> plt.imshow(result, transform=crs, extent=crs.bounds, origin='upper', cmap='RdBu_r')
>>> plt.savefig('tb37v_bmng.png', bbox_inches='tight')
... radius_of_influence=20000, fill_value=None) # doctest: +SKIP
>>> crs = area_def.to_cartopy_crs() # doctest: +SKIP
>>> ax = plt.axes(projection=crs) # doctest: +SKIP
>>> ax.background_img(name='BM') # doctest: +SKIP
>>> plt.imshow(result, transform=crs, extent=crs.bounds, origin='upper', cmap='RdBu_r') # doctest: +SKIP
>>> plt.savefig('tb37v_bmng.png', bbox_inches='tight') # doctest: +SKIP


The above provides you have the Bluemarble background data available in the
Cartopy standard place or in a directory pointed to by the environment
parameter `CARTOPY_USER_BACKGROUNDS`.
Expand Down