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

Allow passing in Red, Green, Blue grids into grdimage #370

Closed
wants to merge 2 commits into from
Closed
Show file tree
Hide file tree
Changes from all 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
37 changes: 31 additions & 6 deletions pygmt/base_plotting.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
Base class with plot generating commands.
Does not define any special non-GMT methods (savefig, show, etc).
"""
import contextlib
import csv
import os
import numpy as np
Expand Down Expand Up @@ -294,27 +295,51 @@ def grdimage(self, grid, **kwargs):
Project grids or images and plot them on maps.

Takes a grid file name or an xarray.DataArray object as input.
Alternatively, pass in a list of red, green, blue grids to be imaged.

Full option list at :gmt-docs:`grdimage.html`

{aliases}

Parameters
----------
grid : str or xarray.DataArray
grid : str, xarray.DataArray or list
The file name of the input grid or the grid loaded as a DataArray.

For plotting RGB grids, pass in a list made up of either file names or
DataArrays to the individual red, green and blue grids.
"""
kwargs = self._preprocess(**kwargs)
kind = data_kind(grid, None, None)

if isinstance(grid, list):
if all([data_kind(g) == "file" for g in grid]):
kind = "file"
grid = " ".join(grid)
elif all([data_kind(g) == "grid" for g in grid]):
kind = "grid"
else:
kind = data_kind(grid)

with Session() as lib:
if kind == "file":
file_context = dummy_context(grid)
file_contexts = [dummy_context(grid)]
elif kind == "grid":
file_context = lib.virtualfile_from_grid(grid)
if isinstance(grid, list):
file_contexts = [
lib.virtualfile_from_grid(grid[0]),
lib.virtualfile_from_grid(grid[1]),
lib.virtualfile_from_grid(grid[2]),
]
else:
file_contexts = [lib.virtualfile_from_grid(grid)]
else:
raise GMTInvalidInput("Unrecognized data type: {}".format(type(grid)))
with file_context as fname:
with contextlib.ExitStack() as stack:
fname = " ".join(
[
stack.enter_context(file_context)
for file_context in file_contexts
]
)
arg_str = " ".join([fname, build_arg_string(kwargs)])
lib.call_module("grdimage", arg_str)

Expand Down
Binary file added pygmt/tests/baseline/test_grdimage_rgb_files.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added pygmt/tests/baseline/test_grdimage_rgb_grid.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
33 changes: 33 additions & 0 deletions pygmt/tests/test_grdimage.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
"""
import numpy as np
import pytest
import xarray as xr

from .. import Figure
from ..exceptions import GMTInvalidInput
Expand Down Expand Up @@ -41,6 +42,38 @@ def test_grdimage_file():
return fig


@pytest.mark.mpl_image_compare
def test_grdimage_rgb_files():
"Plot an image using Red, Green, and Blue file inputs"
fig = Figure()
fig.grdimage(grid=["@earth_relief_60m", "@earth_relief_60m", "@earth_relief_60m"])
Copy link
Member Author

@weiji14 weiji14 Nov 6, 2019

Choose a reason for hiding this comment

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

I don't quite like this unit test example, since they're clearly not RGB grids. Does anyone have a better sample dataset recommendation? Preferably a small one that has a fairly stable URL.

Copy link
Member

Choose a reason for hiding this comment

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

@weiji14 ideally, we should use a dataset that is open-access and published somewhere (figshare, zenodo, etc) or in the public domain. Maybe @djhoese can help us out here 🙂

Copy link

Choose a reason for hiding this comment

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

Hm I can only think of raw data at the moment. What format do you need as input (I'm not super familiar with the limitations of gmt)?

Copy link
Member Author

Choose a reason for hiding this comment

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

GeoTiffs would be nice. Or anything GDAL can read I think. I quite like the zenodo idea, something with a persistent doi. But happy with any other good public domain ones too!

Copy link

Choose a reason for hiding this comment

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

There has to be some landsat data or something. @snowman2 does a lot with rioxarray and may have an idea of a persistent/versioned RGB geotiff. I'll ask the Pytroll core developers too.

Choose a reason for hiding this comment

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

If I understand what you are asking, there is an RGB tif in the example here: https://corteva.github.io/rioxarray/html/examples/COG.html

Copy link
Member Author

Choose a reason for hiding this comment

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

I did look at the rioxarray example a while ago actually, but the test here requires tiffs of individual bands, e.g r.tif, g.tif, b.tif. Landsat on AWS does provide individual bands but they're quite big...

Copy link

Choose a reason for hiding this comment

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

Oh if you need the bands to be in separate tif files then that is going to be harder to track down.

Copy link

@snowman2 snowman2 Nov 25, 2019

Choose a reason for hiding this comment

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

You could just take a regular RGB file and write each band to its own tif file.

import rioxarray
rds = rioxarray,open_rasterio(....)
rdssub = rds.sel(band=1)
rdssub.rio.to_raster(....)

Copy link
Member Author

@weiji14 weiji14 Nov 25, 2019

Choose a reason for hiding this comment

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

Yeah, I thought I might need to create some dummy test files. @seisman just mentioned a similar idea with splitting a raster into R,G,B bands. Just seems like quite a bit of work for a unit test (requires set-up and tear down). What about daily/weather satellites e.g. MODIS/Himawari with a coarse resolution, anyone know where to get a sample?

return fig


@pytest.mark.mpl_image_compare
def test_grdimage_rgb_grid():
"Plot an image using Red, Green, and Blue xarray.DataArray inputs"
red = xr.DataArray(
data=[[128, 0, 0], [128, 0, 0]],
dims=("lat", "lon"),
coords={"lat": [0, 1], "lon": [2, 3, 4]},
)
green = xr.DataArray(
data=[[0, 128, 0], [0, 128, 0]],
dims=("lat", "lon"),
coords={"lat": [0, 1], "lon": [2, 3, 4]},
)
blue = xr.DataArray(
data=[[0, 0, 128], [0, 0, 128]],
dims=("lat", "lon"),
coords={"lat": [0, 1], "lon": [2, 3, 4]},
)

fig = Figure()
fig.grdimage(grid=[red, green, blue], projection="x5c", frame=True)
return fig


def test_grdimage_fails():
"Should fail for unrecognized input"
fig = Figure()
Expand Down