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

seisman · 2023-09-22T05:23:34Z

.tif means a standard TIFF file, while .tiff means a GeoTIFF file.

import pygmt
fig = pygmt.Figure()
fig.coast(region=[0, 10, 0, 10], projection="M10c", frame=True, shorelines="1p")
fig.savefig("map.tiff")

Closes #2658.

Preview: https://pygmt-dev--2698.org.readthedocs.build/en/2698/api/generated/pygmt.Figure.savefig.html

pygmt/figure.py

pygmt/tests/test_figure.py

pygmt/figure.py

pygmt/tests/test_figure.py

Co-authored-by: Yvonne Fröhlich <94163266+yvonnefroehlich@users.noreply.github.com>

pygmt/tests/test_figure.py

weiji14 · 2023-09-24T21:34:47Z

pygmt/tests/test_figure.py

+    geofname = Path("test_figure_savefig_geotiff.tiff")
+    fig.savefig(geofname)
+    assert geofname.exists()
+    assert geofname.with_suffix(".pgw").exists()  # The companion world file


Strange, why is a pgw file generated instead of a tfw file? See https://en.wikipedia.org/wiki/World_file#Filename_extension. A .tif or .tiff file should be linked to .tfw no?

See #2658 (comment) for some tests.

Here, we call gmt psconvert -A -P -W+g without specifying the -T option.

Looking at the psconvert source code (https://github.com/GenericMappingTools/gmt/blob/master/src/psconvert.c#L1008-L1010),

if (!Ctrl->T.active) { /* Set default output device if none is specified */ Ctrl->T.device = (Ctrl->W.warp) ? GS_DEV_PNG : GS_DEV_JPG; /* Lossless PNG if we are making a geotiff in the end */ }

The default format is PNG is -W+g is used. I think GMT fist converts PS to lossless PNG, then write the world file, and calls gdal_translate to combine the PNG and world file into a GeoTIFF file. The extension of the world file is determined from the extension of the raster file (in this case PNG), thus we have a pgw file. This may be an upstream feature or bug.

Actually, since we already have a GeoTIFF file, the world_file is no longer needed. Maybe we should remove it?

Ah ok, thanks for the explanation. You're right that we wouldn't need the pgw file or any world file anyway since the GeoTIFF should already have the metadata embedded inside.

Maybe we should remove it?

Maybe best to add the code to remove the pgw file upstream in GMT when a GeoTIFF is created, or at least rename pgw to tfw? On the PyGMT side, we could match the upstream implementation with GMT <= 6.4.0.

Maybe best to add the code to remove the pgw file upstream in GMT when a GeoTIFF is created, or at least rename pgw to tfw?

Removing pgw makes more sense to me. I'll open an issue report and see what the GMT team thinks about it.

Edit: upstream issue report at GenericMappingTools/gmt#7844

On the PyGMT side, we could match the upstream implementation with GMT <= 6.4.0.

Do you mean we should keep the pgw file? I'm thinking about adding a new parameter worldfile=True to Figure.savefig() which generates a companion world file.

On the PyGMT side, we could match the upstream implementation with GMT <= 6.4.0.

Do you mean we should keep the pgw file? I'm thinking about adding a new parameter worldfile=True to Figure.savefig() which generates a companion world file.

I mean we should match whatever GMT does with the pgw world file later in GMT 6.5.0 (e.g. delete the pgw file), but open up that issue first. The worldfile=True parameter sounds nice actually, this would be similar to GDAL's TFW=YES/NO option when creating GeoTIFFs - https://gdal.org/drivers/raster/gtiff.html#creation-options

Note that GenericMappingTools/gmt#7865 was merged, and the upcoming GMT 6.5.0 should delete the pgw world file if a GeoTIFF output is requested.

pygmt/figure.py

Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com>

seisman · 2023-10-07T03:36:36Z

Ping @GenericMappingTools/pygmt-maintainers for reviews

pygmt/figure.py

weiji14 · 2023-10-08T21:57:49Z

pygmt/tests/test_figure.py

+    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
+
+        # GeoTIFF
+        with rioxarray.open_rasterio(geofname) as xds:
+            assert xds.rio.crs is not None


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:

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.

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.

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:

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.

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.

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.

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):

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.

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.

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

Bottom-right corner

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).

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.

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.

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.

@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.

pygmt/tests/test_figure.py

EJFielding · 2023-10-08T22:42:10Z

In case it helps, here is the GMT6 classic mode script that I used to make the GeoTIFF map (sorry for the ugly csh scripting):

# chose file to plot
set file = $1
set scale = $2
# CSI psvelo output in cm, not m
set GPSscale = 0.2
set Fault = 'Acapulco'
set topoDir = ../../topo
set vecDir = $WORK_DIR/vectors
set GPSdir = $WORK_DIR/modeling/CSI/GPS
set GPSfile = 'Acapulco_v1.psvelo'
set GPSstations = 'Acapulco_v1_sites.txt'

set LLreg = -R-101.0/16.0/-98.5/18.0r
set LLproj = -Ju14/1:2000000

mapproject Acapulco-corn-dd.xy -R-102/-98/10/20 -Ju14/1:1 -F -C > Acapulco-corn-UTM14.xy
# these map units seem to be 100 times greater, probably cm vs. m
set proj = -Jx0.5e-4
set reg = "-R285980/552932/1769965/1990256"
set projkm = -Jx0.5e-1
set regkm = "-R285.980/552.932/1769.965/1990.256"

# new map view with vectors and error ellipses
set mapOut = $Fault-slip-map5.eps

set slipd = ${Fault}.slipdir
set fault = $Fault.fault

set out = $file-Acapulco-map5.eps
set in_cpt = polar
set cpt = co_unw-red_blue-m.cpt
set col_scale = "-T-$scale/$scale/0.01"

set resol = ""

# set for tight bounding box
gmt set PS_MEDIA letter
gmt set FORMAT_GEO_OUT Ddd.x

gmt makecpt -C$in_cpt $col_scale -Z -D > $cpt
gmt psbasemap $LLproj $LLreg  -Bxy0.2ag --MAP_FRAME_TYPE=inside -K -V -P > $out

# fault slip rectangles colored by their -Z values
gmt psxy $file.dat $LLproj $LLreg -C$cpt -L -W0.01c -O -K >> $out

# better to have no fill for GeoTIFF
gmt pscoast $LLproj $LLreg -N1/0/0/0 -I1/0/0/255 -Wthick,black -Df -O -K >> $out


# GPS stations
gmt psxy $GPSdir/$GPSstations $LLproj $LLreg -fi2x,1y -h1 -St0.4 -Wthick,yellow -Gyellow -V -O -K >> $out
gmt psxy $GPSdir/$GPSstations $LLproj $LLreg -fi2x,1y -h1 -St0.3 -Wthin,black -Gblack -V -O -K >> $out

gmt psvelo GPS_data_transformed_Acapulco.psvelo $LLproj $LLreg -Se$GPSscale/0.95 -A0.07/0.20/0.15 -G0/0/0  -P -K -O -V >> $out
gmt psvelo GPS_synthetic_Acapulco.psvelo $LLproj $LLreg -Se$GPSscale/0.95+f0p -A0.07/0.20/0.15 -Ggreen  -P -K -O -V >> $out

# fault rectangle on top--model fault already in UTM km
gmt psxy $vecDir/M7_1_SSN_MX_hypocenter_early_verificado.llz $LLproj $LLreg -Sa0.45 -W2/0/0/0 -Ggreen -O -K >> $out

gmt psscale -Dx0.2c/1.5c+w2.5c/0.4+e -C$cpt -Bxaf+1"Slip (m)" -O -V  >> $out

gv $out

gmt psconvert -A -P -Tt -W+g $out

Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com>

weiji14

An extra check for the Affine transformation, otherwise ok.

pygmt/tests/test_figure.py

Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com>

EJFielding · 2023-12-29T18:41:32Z

Thanks to @seisman and @weiji14 for adding this GeoTIFF function to PyGMT! I am now using it all the time, and it is working great.

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

41767a2

seisman added the enhancement Improving an existing feature label Sep 22, 2023

seisman added this to the 0.11.0 milestone Sep 22, 2023

seisman mentioned this pull request Sep 22, 2023

Figure.psconvert: Add support for GeoTIFF output #2658

Closed

seisman added 6 commits September 22, 2023 13:35

Test the .tiff extension

cd9918f

GeoTIFF doesn't need the -T option

f2348cc

Add a new test to check if .tiff is a GeoTIFF file

921912d

Revert the changes in test_figure_savefig_exists

a83862f

Simplify the test

51c8dbe

Fix a linting issue

cf771cc

seisman marked this pull request as ready for review September 22, 2023 15:32

seisman added needs review This PR has higher priority and needs review. feature Brand new feature and removed enhancement Improving an existing feature labels Sep 22, 2023

seisman requested a review from a team September 24, 2023 01:54

yvonnefroehlich reviewed Sep 24, 2023

View reviewed changes

pygmt/figure.py Outdated Show resolved Hide resolved

pygmt/tests/test_figure.py Outdated Show resolved Hide resolved

pygmt/figure.py Outdated Show resolved Hide resolved

pygmt/tests/test_figure.py Outdated Show resolved Hide resolved

Fix typos

3303ba4

Co-authored-by: Yvonne Fröhlich <94163266+yvonnefroehlich@users.noreply.github.com>

weiji14 reviewed Sep 24, 2023

View reviewed changes

Apply suggestions from code review

4fd2f50

Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com>

seisman mentioned this pull request Sep 25, 2023

Incorrect world file name when generating GeoTIFF file GenericMappingTools/gmt#7844

Closed

seisman added 7 commits September 29, 2023 17:36

Merge branch 'main' into geotiff

9e77ef3

Merge branch 'main' into geotiff

d5e6617

Merge remote-tracking branch 'origin/geotiff' into geotiff

f1f41bf

Delete the .pgw if exists

252514b

Merge branch 'main' into geotiff

c4d480e

Merge branch 'main' into geotiff

2137682

Merge branch 'main' into geotiff

d57d8e0

weiji14 reviewed Oct 8, 2023

View reviewed changes

pygmt/tests/test_figure.py Show resolved Hide resolved

weiji14 reviewed Oct 8, 2023

View reviewed changes

pygmt/tests/test_figure.py Show resolved Hide resolved

seisman and others added 7 commits October 9, 2023 10:10

Improve the comment about removing .pgw file

44bfcfe

Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com>

Check bounds in the test

e7c2b36

Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com>

Merge branch 'main' into geotiff

9127756

Fix typos and linting issues

26837e8

Merge branch 'main' into geotiff

d5ececc

Merge branch 'main' into geotiff

05ef422

Fix the bounds and shape for geotiff files

71611df

seisman requested a review from a team October 18, 2023 09:04

Merge branch 'main' into geotiff

916b4a8

weiji14 approved these changes Oct 24, 2023

View reviewed changes

pygmt/tests/test_figure.py Show resolved Hide resolved

pygmt/tests/test_figure.py Show resolved Hide resolved

pygmt/tests/test_figure.py Show resolved Hide resolved

seisman and others added 2 commits October 24, 2023 21:13

Check Affine transformation

6a1caa3

Co-authored-by: Wei Ji <23487320+weiji14@users.noreply.github.com>

Merge branch 'main' into geotiff

da9ca29

seisman removed the needs review This PR has higher priority and needs review. label Oct 24, 2023

seisman merged commit 904553d into main Oct 24, 2023
16 checks passed

seisman deleted the geotiff branch October 24, 2023 13:44

seisman mentioned this pull request Oct 24, 2023

Figure.savefig: Add the 'worldfile' parameter to write a companion world file for raster images #2766

Merged

7 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

seisman commented Sep 22, 2023 •

edited by yvonnefroehlich

Loading

weiji14 Sep 24, 2023

seisman Sep 25, 2023

weiji14 Sep 25, 2023

seisman Sep 25, 2023 •

edited

Loading

weiji14 Sep 25, 2023 •

edited

Loading

weiji14 Oct 8, 2023

seisman commented Oct 7, 2023

weiji14 Oct 8, 2023

weiji14 Oct 8, 2023

seisman Oct 9, 2023

weiji14 Oct 9, 2023 •

edited

Loading

seisman Oct 9, 2023

weiji14 Oct 9, 2023

seisman Oct 9, 2023

EJFielding Oct 9, 2023

EJFielding Oct 9, 2023 •

edited

Loading

seisman Oct 18, 2023

EJFielding commented Oct 8, 2023

weiji14 left a comment

EJFielding commented Dec 29, 2023

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

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

Conversation

seisman commented Sep 22, 2023 • edited by yvonnefroehlich Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

seisman Sep 25, 2023 • edited Loading

Choose a reason for hiding this comment

weiji14 Sep 25, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

seisman commented Oct 7, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

weiji14 Oct 9, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

EJFielding Oct 9, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

EJFielding commented Oct 8, 2023

weiji14 left a comment

Choose a reason for hiding this comment

EJFielding commented Dec 29, 2023

seisman commented Sep 22, 2023 •

edited by yvonnefroehlich

Loading

seisman Sep 25, 2023 •

edited

Loading

weiji14 Sep 25, 2023 •

edited

Loading

weiji14 Oct 9, 2023 •

edited

Loading

EJFielding Oct 9, 2023 •

edited

Loading