From 9d3831632c2a0fb7ba10f6cb2ade87343128ab17 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Fri, 25 Aug 2023 18:36:33 +1200 Subject: [PATCH 1/9] Add gallery example for plotting an RGB image from an xarray.DataArray MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Gallery example using pygmt.Figure.grdimage to plot an RGB image from a 3-band GeoTIFF loaded into an xarray.DataArray via rioxarray.open_rasterio. Example is over Lāhainā, Hawai'i on 9 Aug 2023. --- examples/gallery/images/rgb_image.py | 37 ++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) create mode 100644 examples/gallery/images/rgb_image.py diff --git a/examples/gallery/images/rgb_image.py b/examples/gallery/images/rgb_image.py new file mode 100644 index 00000000000..61040eb7403 --- /dev/null +++ b/examples/gallery/images/rgb_image.py @@ -0,0 +1,37 @@ +""" +RGB Image +--------- +The :meth:`pygmt.Figure.grdimage` method can be used to plot Red, Green, Blue +(RGB) images, or any 3-band false color combination. Here, we'll use +:meth:`rioxarray.open_rasterio` to read a GeoTIFF file into an +:class:`xarray.DataArray` format, and plot it on a map. + +The example below shows a Worldview 2 satellite image over +`Lāhainā, Hawai'i during the August 2023 wildfires +`_. +Data is sourced from a Cloud-Optimized GeoTIFF file hosted on +`OpenAerialMap `_ under a +`CC BY-NC 4.0 `_ license. +""" +import pygmt +import rioxarray + +############################################################################### +# Read 3-band data from GeoTIFF into an xarray.DataArray object +image = rioxarray.open_rasterio( + filename="https://oin-hotosm.s3.us-east-1.amazonaws.com/64d6a49a19cb3a000147a65b/0/64d6a49a19cb3a000147a65c.tif", + overview_level=5, +) +# Subset to area of Lāhainā in EPSG:32604 coordinates +image = image.rio.clip_box(minx=738000, maxx=755000, miny=2300000, maxy=2318000) +image + +############################################################################### +# Plot the RGB imagery +fig = pygmt.Figure() +fig.grdimage( + grid=image, + projection="x1:100000", + frame=[r"WSne+tL@!a\225hain@!a\225, Hawai\047i on 9 Aug 2023", "af"], +) +fig.show() From 60d1be38303a23a1d239076276da82f4a7873400 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Sat, 26 Aug 2023 13:15:11 +1200 Subject: [PATCH 2/9] Add abbreviation for Cloud-Optimized GeoTIFF Co-Authored-By: Michael Grund <23025878+michaelgrund@users.noreply.github.com> --- examples/gallery/images/rgb_image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/gallery/images/rgb_image.py b/examples/gallery/images/rgb_image.py index 61040eb7403..6929c880a8c 100644 --- a/examples/gallery/images/rgb_image.py +++ b/examples/gallery/images/rgb_image.py @@ -9,7 +9,7 @@ The example below shows a Worldview 2 satellite image over `Lāhainā, Hawai'i during the August 2023 wildfires `_. -Data is sourced from a Cloud-Optimized GeoTIFF file hosted on +Data is sourced from a Cloud-Optimized GeoTIFF (COG) file hosted on `OpenAerialMap `_ under a `CC BY-NC 4.0 `_ license. """ From b516c194da89836ecc83b745775ffe41bfb96b88 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Sat, 26 Aug 2023 13:21:02 +1200 Subject: [PATCH 3/9] Set title font to Times-Roman MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Makes the unequal spacing before and after the ā less obvious. --- examples/gallery/images/rgb_image.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/examples/gallery/images/rgb_image.py b/examples/gallery/images/rgb_image.py index 6929c880a8c..64101c8b257 100644 --- a/examples/gallery/images/rgb_image.py +++ b/examples/gallery/images/rgb_image.py @@ -29,9 +29,11 @@ ############################################################################### # Plot the RGB imagery fig = pygmt.Figure() -fig.grdimage( - grid=image, - projection="x1:100000", - frame=[r"WSne+tL@!a\225hain@!a\225, Hawai\047i on 9 Aug 2023", "af"], -) +with pygmt.config(FONT_TITLE="Times-Roman"): # Set title font to Times-Roman + fig.grdimage( + grid=image, + projection="x1:100000", + frame=[r"WSne+tL@!a\225hain@!a\225, Hawai\047i on 9 Aug 2023", "af"], + ) +fig.savefig("Lāhainā.png") fig.show() From e7b66e6cb4c83b8854778addab4ba0ae13ff492b Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Sat, 26 Aug 2023 14:10:18 +1200 Subject: [PATCH 4/9] =?UTF-8?q?Change=20apostrophe=20to=20Okina=20symbol?= =?UTF-8?q?=20=CA=BB?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit See https://en.wikipedia.org/wiki/%CA%BBOkina. Using octal code 140 instead of 047 on the plot title for Hawaiʻi, so that it looks like a 6 instead of a 9, but unsure if this is still the correct Okina symbol. --- examples/gallery/images/rgb_image.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/gallery/images/rgb_image.py b/examples/gallery/images/rgb_image.py index 64101c8b257..2d1f2a54503 100644 --- a/examples/gallery/images/rgb_image.py +++ b/examples/gallery/images/rgb_image.py @@ -7,7 +7,7 @@ :class:`xarray.DataArray` format, and plot it on a map. The example below shows a Worldview 2 satellite image over -`Lāhainā, Hawai'i during the August 2023 wildfires +`Lāhainā, Hawaiʻi during the August 2023 wildfires `_. Data is sourced from a Cloud-Optimized GeoTIFF (COG) file hosted on `OpenAerialMap `_ under a @@ -33,7 +33,7 @@ fig.grdimage( grid=image, projection="x1:100000", - frame=[r"WSne+tL@!a\225hain@!a\225, Hawai\047i on 9 Aug 2023", "af"], + frame=[r"WSne+tL@!a\225hain@!a\225, Hawai\140i on 9 Aug 2023", "af"], ) fig.savefig("Lāhainā.png") fig.show() From 9085f6586140a616d81ded1ff4050bed6377c805 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Sat, 26 Aug 2023 15:01:33 +1200 Subject: [PATCH 5/9] Inline comment on map scale projection The 1:100000 scale means 1 centimetre on the map is equivalent to 1 kilometre on the ground. --- examples/gallery/images/rgb_image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/gallery/images/rgb_image.py b/examples/gallery/images/rgb_image.py index 2d1f2a54503..eea1e729502 100644 --- a/examples/gallery/images/rgb_image.py +++ b/examples/gallery/images/rgb_image.py @@ -32,7 +32,7 @@ with pygmt.config(FONT_TITLE="Times-Roman"): # Set title font to Times-Roman fig.grdimage( grid=image, - projection="x1:100000", + projection="x1:100000", # map scale where 1 cm on map equals 1 km on the ground frame=[r"WSne+tL@!a\225hain@!a\225, Hawai\140i on 9 Aug 2023", "af"], ) fig.savefig("Lāhainā.png") From 4fd90bee2a1bf306d22d31a24054042050ac6eed Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Sat, 26 Aug 2023 20:11:58 +1200 Subject: [PATCH 6/9] =?UTF-8?q?Use=20=C2=AF=20and=20`=20instead=20of=20\22?= =?UTF-8?q?5=20and=20\140=20The=20non-octal=20code=20versions=20are=20much?= =?UTF-8?q?=20easier=20to=20see.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Dongdong Tian --- examples/gallery/images/rgb_image.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/gallery/images/rgb_image.py b/examples/gallery/images/rgb_image.py index eea1e729502..631155975fc 100644 --- a/examples/gallery/images/rgb_image.py +++ b/examples/gallery/images/rgb_image.py @@ -33,7 +33,6 @@ fig.grdimage( grid=image, projection="x1:100000", # map scale where 1 cm on map equals 1 km on the ground - frame=[r"WSne+tL@!a\225hain@!a\225, Hawai\140i on 9 Aug 2023", "af"], + frame=[r"WSne+tL@!a¯hain@!a¯, Hawai`i on 9 Aug 2023", "af"], ) -fig.savefig("Lāhainā.png") fig.show() From 3c0712d5b9b1e02e1572d21dcd4ff3978f3b7324 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Wed, 30 Aug 2023 15:38:24 +1200 Subject: [PATCH 7/9] Force loading dataarray into memory Opening the GeoTIFF in a with statement, and loading the RGB image data fully into memory using `image.load()` as suggested in https://github.com/corteva/rioxarray/issues/550#issuecomment-1216690984, which should fix the segfault hopefully. --- examples/gallery/images/rgb_image.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/gallery/images/rgb_image.py b/examples/gallery/images/rgb_image.py index 631155975fc..78df1c3e1e7 100644 --- a/examples/gallery/images/rgb_image.py +++ b/examples/gallery/images/rgb_image.py @@ -18,12 +18,13 @@ ############################################################################### # Read 3-band data from GeoTIFF into an xarray.DataArray object -image = rioxarray.open_rasterio( +with rioxarray.open_rasterio( filename="https://oin-hotosm.s3.us-east-1.amazonaws.com/64d6a49a19cb3a000147a65b/0/64d6a49a19cb3a000147a65c.tif", overview_level=5, -) -# Subset to area of Lāhainā in EPSG:32604 coordinates -image = image.rio.clip_box(minx=738000, maxx=755000, miny=2300000, maxy=2318000) +) as img: + # Subset to area of Lāhainā in EPSG:32604 coordinates + image = img.rio.clip_box(minx=738000, maxx=755000, miny=2300000, maxy=2318000) + image = image.load() # force loading dataarray into memory image ############################################################################### From 00b427b4d83240922edeab369fb1404d03309445 Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Thu, 31 Aug 2023 09:55:12 +1200 Subject: [PATCH 8/9] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Yvonne Fröhlich <94163266+yvonnefroehlich@users.noreply.github.com> --- examples/gallery/images/rgb_image.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/examples/gallery/images/rgb_image.py b/examples/gallery/images/rgb_image.py index 78df1c3e1e7..47af155b2cf 100644 --- a/examples/gallery/images/rgb_image.py +++ b/examples/gallery/images/rgb_image.py @@ -3,7 +3,7 @@ --------- The :meth:`pygmt.Figure.grdimage` method can be used to plot Red, Green, Blue (RGB) images, or any 3-band false color combination. Here, we'll use -:meth:`rioxarray.open_rasterio` to read a GeoTIFF file into an +:py:func:`rioxarray.open_rasterio` to read a GeoTIFF file into an :class:`xarray.DataArray` format, and plot it on a map. The example below shows a Worldview 2 satellite image over @@ -17,23 +17,24 @@ import rioxarray ############################################################################### -# Read 3-band data from GeoTIFF into an xarray.DataArray object +# Read 3-band data from GeoTIFF into an xarray.DataArray object: with rioxarray.open_rasterio( filename="https://oin-hotosm.s3.us-east-1.amazonaws.com/64d6a49a19cb3a000147a65b/0/64d6a49a19cb3a000147a65c.tif", overview_level=5, ) as img: # Subset to area of Lāhainā in EPSG:32604 coordinates image = img.rio.clip_box(minx=738000, maxx=755000, miny=2300000, maxy=2318000) - image = image.load() # force loading dataarray into memory + image = image.load() # Force loading the DataArray into memory image ############################################################################### -# Plot the RGB imagery +# Plot the RGB imagery: fig = pygmt.Figure() with pygmt.config(FONT_TITLE="Times-Roman"): # Set title font to Times-Roman fig.grdimage( grid=image, - projection="x1:100000", # map scale where 1 cm on map equals 1 km on the ground + # Use a map scale where 1 cm on the map equals 1 km on the ground + projection="x1:100000", frame=[r"WSne+tL@!a¯hain@!a¯, Hawai`i on 9 Aug 2023", "af"], ) fig.show() From 2553ccee6be3caf3864d834291627ddeb8a9ec3a Mon Sep 17 00:00:00 2001 From: Wei Ji <23487320+weiji14@users.noreply.github.com> Date: Thu, 31 Aug 2023 09:57:00 +1200 Subject: [PATCH 9/9] Lint --- examples/gallery/images/rgb_image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/gallery/images/rgb_image.py b/examples/gallery/images/rgb_image.py index 47af155b2cf..9d9506b2ffb 100644 --- a/examples/gallery/images/rgb_image.py +++ b/examples/gallery/images/rgb_image.py @@ -34,7 +34,7 @@ fig.grdimage( grid=image, # Use a map scale where 1 cm on the map equals 1 km on the ground - projection="x1:100000", + projection="x1:100000", frame=[r"WSne+tL@!a¯hain@!a¯, Hawai`i on 9 Aug 2023", "af"], ) fig.show()