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

Figure.savefig: Support generating GeoTIFF file (with extension '.tiff') #2698

Merged
merged 26 commits into from
Oct 24, 2023
Merged
Show file tree
Hide file tree
Changes from 24 commits
Commits
Show all changes
26 commits
Select commit Hold shift + click to select a range
41767a2
Figure.savefig: Support generating GeoTIFF file (with extension '.tiff')
seisman Sep 22, 2023
cd9918f
Test the .tiff extension
seisman Sep 22, 2023
f2348cc
GeoTIFF doesn't need the -T option
seisman Sep 22, 2023
921912d
Add a new test to check if .tiff is a GeoTIFF file
seisman Sep 22, 2023
a83862f
Revert the changes in test_figure_savefig_exists
seisman Sep 22, 2023
51c8dbe
Simplify the test
seisman Sep 22, 2023
cf771cc
Fix a linting issue
seisman Sep 22, 2023
3303ba4
Fix typos
seisman Sep 24, 2023
4fd2f50
Apply suggestions from code review
seisman Sep 24, 2023
9e77ef3
Merge branch 'main' into geotiff
seisman Sep 29, 2023
d5e6617
Merge branch 'main' into geotiff
seisman Oct 4, 2023
f1f41bf
Merge remote-tracking branch 'origin/geotiff' into geotiff
seisman Oct 4, 2023
252514b
Delete the .pgw if exists
seisman Oct 4, 2023
c4d480e
Merge branch 'main' into geotiff
seisman Oct 4, 2023
2137682
Merge branch 'main' into geotiff
seisman Oct 5, 2023
d57d8e0
Merge branch 'main' into geotiff
seisman Oct 7, 2023
44bfcfe
Improve the comment about removing .pgw file
seisman Oct 9, 2023
e7c2b36
Check bounds in the test
seisman Oct 9, 2023
9127756
Merge branch 'main' into geotiff
seisman Oct 9, 2023
26837e8
Fix typos and linting issues
seisman Oct 9, 2023
d5ececc
Merge branch 'main' into geotiff
seisman Oct 12, 2023
05ef422
Merge branch 'main' into geotiff
seisman Oct 18, 2023
71611df
Fix the bounds and shape for geotiff files
seisman Oct 18, 2023
916b4a8
Merge branch 'main' into geotiff
seisman Oct 18, 2023
6a1caa3
Check Affine transformation
seisman Oct 24, 2023
da9ca29
Merge branch 'main' into geotiff
seisman Oct 24, 2023
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
41 changes: 28 additions & 13 deletions pygmt/figure.py
Original file line number Diff line number Diff line change
Expand Up @@ -257,12 +257,18 @@ def savefig(
"""
Save the figure to a file.

This method implements a matplotlib-like interface for
:meth:`pygmt.Figure.psconvert`.
Supported file formats and their extensions:

Supported formats: PNG (``.png``), JPEG (``.jpg`` or ``.jpeg``),
PDF (``.pdf``), BMP (``.bmp``), TIFF (``.tif``), EPS (``.eps``), and
KML (``.kml``). The KML output generates a companion PNG file.
- PNG (``.png``)
- JPEG (``.jpg`` or ``.jpeg``)
- PDF (``.pdf``)
- BMP (``.bmp``)
- TIFF (``.tif``)
- GeoTIFF (``.tiff``)
- EPS (``.eps``)
- KML (``.kml``)

For KML format, a companion PNG file is also generated.

You can pass in any keyword arguments that
:meth:`pygmt.Figure.psconvert` accepts.
Expand All @@ -279,10 +285,10 @@ def savefig(
If ``True``, will crop the figure canvas (page) to the plot area.
anti_alias: bool
If ``True``, will use anti-aliasing when creating raster images
(PNG, JPG, TIFF). More specifically, it passes arguments ``t2``
and ``g2`` to the ``anti_aliasing`` parameter of
:meth:`pygmt.Figure.psconvert`. Ignored if creating vector
graphics.
(BMP, PNG, JPEG, TIFF, and GeoTIFF). More specifically, it passes
the arguments ``"t2"`` and ``"g2"`` to the ``anti_aliasing``
parameter of :meth:`pygmt.Figure.psconvert`. Ignored if creating
vector graphics.
show: bool
If ``True``, will open the figure in an external viewer.
dpi : int
Expand All @@ -301,15 +307,20 @@ def savefig(
"bmp": "b",
"eps": "e",
"tif": "t",
"tiff": None, # GeoTIFF doesn't need the -T option
"kml": "g",
}

fname = Path(fname)
prefix, suffix = fname.with_suffix("").as_posix(), fname.suffix
ext = suffix[1:].lower() # Remove the . and normalize to lowercase
# alias jpeg to jpg
if ext == "jpeg":

if ext == "jpeg": # Alias jpeg to jpg
ext = "jpg"
elif ext == "tiff": # GeoTIFF
kwargs["W"] = "+g"
elif ext == "kml": # KML
kwargs["W"] = "+k"

if ext not in fmts:
if ext == "ps":
Expand All @@ -328,11 +339,15 @@ def savefig(
if anti_alias:
kwargs["Qt"] = 2
kwargs["Qg"] = 2
if ext == "kml":
kwargs["W"] = "+k"

self.psconvert(prefix=prefix, fmt=fmt, crop=crop, **kwargs)

# Remove the .pgw world file if exists
# Not necessary after GMT 6.5.0.
# See upstream fix https://github.com/GenericMappingTools/gmt/pull/7865
if ext == "tiff" and fname.with_suffix(".pgw").exists():
fname.with_suffix(".pgw").unlink()

# Rename if file extension doesn't match the input file suffix
if ext != suffix[1:]:
fname.with_suffix("." + ext).rename(fname)
Expand Down
54 changes: 54 additions & 0 deletions pygmt/tests/test_figure.py
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,60 @@ def test_figure_savefig_exists():
fname.unlink()


def test_figure_savefig_geotiff():
"""
Make sure .tif generates a normal TIFF file and .tiff generates a GeoTIFF
file.
"""
fig = Figure()
fig.basemap(region=[0, 10, 0, 10], projection="M10c", frame=True)

# Save as GeoTIFF
geofname = Path("test_figure_savefig_geotiff.tiff")
fig.savefig(geofname)
assert geofname.exists()
# The .pgw should not exist
assert not geofname.with_suffix(".pgw").exists()

# Save as TIFF
fname = Path("test_figure_savefig_tiff.tif")
fig.savefig(fname)
assert fname.exists()

# Check if a TIFF is georeferenced or not
try:
# pylint: disable=import-outside-toplevel
import rioxarray
from rasterio.errors import NotGeoreferencedWarning
seisman marked this conversation as resolved.
Show resolved Hide resolved

# GeoTIFF
with rioxarray.open_rasterio(geofname) as xds:
assert xds.rio.crs is not None
Copy link
Member

Choose a reason for hiding this comment

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

The output CRS from this 10cm Mercator projection map is:

PROJCS["unknown",GEOGCS["unknown",DATUM["unknown",SPHEROID["unknown",6378137,298.257220143428]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",5],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

I tried dragging and dropping the test_figure_savefig_geotiff.tiff into QGIS, and the projection showed up as 'unknown' (same as rioxarray above), but it defaulted to plotting on longitude 0-10, latitude 0-10:

image

However, the blue gridlines (actual longitude/latitude lines from QGIS) don't match up with the GeoTIFF (from psconvert), see the top left corner around 10°N, where the black/white zebra grid markings don't line up with the blue lines.

image

Not sure what the default projection is that is output from psconvert, but I think we should make sure there isn't a bug somewhere. Maybe we should try a few other projections besides M10c and see if the output is ok?

@EJFielding, could you share with us the projection (-J) parameter you used to generate the overlay plot over Mexico at #2658 (comment)? Just wanted to know how you're aligning the grids properly.

Copy link
Member

Choose a reason for hiding this comment

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

Oh wait nevermind, I think I had a cached tiff.aux.xml file from trying a few other projections like -JEPSG:3857, and it was rendering the GeoTIFF in the wrong projection. Deleting all the files and resaving them seems to work better now:

image

There's still a 900m gap between the 9°N blue longitude line and the black/white zebra grid border, but I think that's just because of the default dpi resolution being a bit coarse.

Copy link
Member Author

Choose a reason for hiding this comment

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

Does the frame make trouble here? The psconvert -W (https://docs.generic-mapping-tools.org/latest/psconvert.html) says:

Be aware, however, that different results are obtained depending on the image contents and if the -B option has been used or not. The trouble with the -B option is that it creates a frame and very likely its annotations. That introduces pixels outside the map data extent, and therefore the map extents estimation will be wrong. To avoid this problem use --MAP_FRAME_TYPE=inside option which plots all annotations and ticks inside the image and therefore does not compromise the coordinate computations. Pay attention also to the cases when the plot has any of the sides with whites only because than the algorithm will fail miserably as those whites will be eaten by the Ghostscript. In that case you really must use -B or use a slightly off-white color.

Copy link
Member

@weiji14 weiji14 Oct 9, 2023

Choose a reason for hiding this comment

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

Assuming that the -B means -B from grdimage/basemap/etc and not psconvert... I tried:

with pygmt.config(MAP_FRAME_TYPE="inside"):
        fig.basemap(region=[0, 10, 0, 10], projection="M10c", frame="afg")

And it does look a bit better on the latitude (y-axis) side on the top left corner. The blue line and black lines align well on the latitude grid.

image

But on the longitude (x-axis) side on the bottom right, there still seems to be some offset (~900m) with the blue line not aligning with the black line (and actually, there's some offset from the 0deg latitude origin too):

image

Given that the world file's georeferencing starts from the top-left corner (see https://en.wikipedia.org/wiki/World_file), it kinda makes sense that the top left corner is more aligned, and the offset/distortion is greater at the bottom right corner. Not sure if there's much we can do about this though.

Copy link
Member Author

Choose a reason for hiding this comment

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

Could you please try to produce a figure using Figure.coast or Figure.grdimage without the frame parameter? The frame still has a thickness, which may affect the boundingbox or something.

Copy link
Member

Choose a reason for hiding this comment

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

Still having an offset without the frame:

fig = Figure()
with pygmt.config(MAP_FRAME_TYPE="inside"):
    fig.coast(shorelines=True, region=[0, 10, 0, 10], projection="M10c")
# Save as GeoTIFF
geofname = Path("test_figure_savefig_geotiff.tiff")
fig.savefig(geofname)

Top-left corner

image

Bottom-right corner

image

I measured the size of one of the gray pixels and it was ~900+m, same as the offset. So it might be a case of gridline vs pixel registration (xref #487)? Though I'd think the offset would be half a pixel instead (~450m).

Copy link
Member Author

Choose a reason for hiding this comment

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

So it might be a case of gridline vs pixel registration (xref #487)?

But this is not a grid, so can't be the gridline/pixel registration issue. Anyway, if something doesn't make sense, then it must be an upstream bug.

Choose a reason for hiding this comment

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

I found that the GeoTIFF images I was making with GMT classic looked better with the inside frame. I don't remember if having the frame outside made a difference in the geolocation.

I am not sure what the conventions are for "world" files, but I know that GDAL has the reference location for the image as the top-left corner of the top-left pixel, unless you change the registration. This is the same as "pixel" registration in GMT.

Copy link

@EJFielding EJFielding Oct 9, 2023

Choose a reason for hiding this comment

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

The other thing to realize for checking image geolocation is that the coastline database used in GMT is something like 40 years old, and not very accurate. I would not rely on the GMT coastline to check the geolocations. In some places the GMT coastline is several hundred meters different from Google Maps or OpenStreetMap.

Copy link
Member Author

Choose a reason for hiding this comment

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

@weiji14 When you have time, it would be better to submit an upstream issue for more in-depth discussions about the possible misalignments. I believe we can do nothing on the PyGMT side.

seisman marked this conversation as resolved.
Show resolved Hide resolved
npt.assert_allclose(
actual=xds.rio.bounds(),
desired=(
-661136.0621116752,
-54631.82709660966,
592385.4459661598,
1129371.7360144067,
),
)
assert xds.rio.shape == (1257, 1331)
seisman marked this conversation as resolved.
Show resolved Hide resolved
# TIFF
with pytest.warns(expected_warning=NotGeoreferencedWarning) as record:
with rioxarray.open_rasterio(fname) as xds:
assert xds.rio.crs is None
seisman marked this conversation as resolved.
Show resolved Hide resolved
npt.assert_allclose(
actual=xds.rio.bounds(), desired=(0.0, 0.0, 1331.0, 1257.0)
)
assert xds.rio.shape == (1257, 1331)
seisman marked this conversation as resolved.
Show resolved Hide resolved
assert len(record) == 1
except ImportError:
pass
geofname.unlink()
fname.unlink()


def test_figure_savefig_directory_nonexists():
"""
Make sure that Figure.savefig() raises a FileNotFoundError when the parent
Expand Down