Is it possible to use WMS or WMTS with PyGMT, i.e. to add layers to a map similarly to Figure.tilemap() for XYZ tiles?
It seems GMT.jl supports WMS at least (Web Map Service), but I’m not able to find any information for PyGMT
Cheers!
PyGMT has a tiles module too. Search the examples or tutorial pages. Not sure where.
Thanks. I’m however not able to find any support for WMS/WMTS services, only XYZ tiles as described here pygmt.Figure.tilemap — PyGMT and Tile maps — PyGMT . The contextily package, which can be used as a source, is only supporting XYZ tiles. Or am I missing something obvious?
No. I think you are right.
seems like not a big deal to retrieve map images using gdal and plot it using pygmt. not nearly as nice as GMT.jl, but still does the job:
import pygmt
from osgeo import gdal
gdal.UseExceptions() # to suppress an annoying warning message
tiff = 'SENTINEL_image.tiff'
wms = "WMS:http://tiles.maps.eox.at/wms?"
ds = gdal.Open( wms )
layers = ds.GetMetadata( "SUBDATASETS" )
#layers.keys() to check layer names
#layers.values() to check layer URLs
layer = layers['SUBDATASET_3_NAME']
ds = None # close
ds = gdal.Open( layer )
Res = 1500 # meters
Res = Res / 111195 # 111195 m/degree for spherical Earth
projWin = [-10, 44, -5, 37] # [ ulx, uly, lrx, lry ]
gdal.Translate( tiff, ds, projWin=projWin, xRes=Res, yRes=Res )
ds = None # close
fig = pygmt.Figure()
fig.grdimage( tiff, frame = ['WSen', 'xa', 'ya'] )
fig.savefig( "SENTINEL_image.png", show=True )
import pygmt
from osgeo import gdal
gdal.UseExceptions() # to suppress an annoying warning message
tiff = 'MODIS_image.tiff'
wms = "WMS:https://gibs-c.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi"
ds = gdal.Open( wms )
layers = ds.GetMetadata( "SUBDATASETS" )
for value in layers.values():
if 'MODIS_Terra_CorrectedReflectance_TrueColor' in value:
break
time="2021-10-29"
layer = f"{value}&TIME={time}"
ds = None # close
ds = gdal.Open( layer )
Res = 1500 # meters
Res = Res / 111195 # 111195 m/degree for spherical Earth
projWin = [9, 43, 22, 32] # [ ulx, uly, lrx, lry ]
gdal.Translate( tiff, ds, projWin=projWin, xRes=Res, yRes=Res )
ds = None # close
fig = pygmt.Figure()
fig.grdimage( tiff, frame = ['WSen', 'xa', 'ya'] )
fig.savefig("MODIS_image.png", show=True)
That’s basically what the GMT.jl’s wms functions do under the hood.
Thanks, this is great! I mark this as the solution.
It would still be nice to have higher-level functionality for WMS, so I created a feature request (Feature request: add support for WMS and WMTS · Issue #3167 · GenericMappingTools/pygmt · GitHub). I’m also wondering if it’s possible to do this without using gdal (e.g. owslib, rasterio), or at least keep everything in memory and avoid having to write the tiff to disk. I’ll post here if I find a solution when I have time to look deeper into that.
personally I think the solution using gdal is good enough and (py)gmt developers can use their time and skills for fixing many more important issues.
storing tiff image locally enables reuse. This appears absolutely essential when e.g. doing multiple adjustments of the same map, as the locally stored image can be reused to avoid unnecessary downloading of the same image data from the provider over and over again.
Reusing the tiff also allows your plotting script to run faster.
… you can call the GMT.jl wms functions from Python. They do not need anything else because GDAL is already a GMT dependency and all is done in memory.
People also seem thinking the online data sources are with us forever or neither the user nor the resources ever go offline.
I’ve witnessed occasions of online resources disappearing without traces or becoming classified. Storing your copy of the (image) data may become important before you know it.
I am far from thinking this is absolutely needed for every possible project, but storing files locally solves this “by default”.
how does this work, are there any examples?
does mixing pygmt and GMT.lj calls work for plotting? i.e. will GMT.lj calls plot into the existing pygmt Figure() (i suspect not but asking never harms)
thanks
There is this juliacall package that lets call Julia from Python, and I have shown an example in that tilemap topic but I have not much more experience using it.
PyGMT and GMT.jl plots cannot be mixed. While GMT.jl allows modern mode it tries to avoid it when possible. On the other hand PyGMT only allows modern mode. Each modern mode session has its own session dir so no communication is possible berween two sessions.
thanks!
It turned out super convenient to use pygmt.clib package for creating GMT session using clib.Session() and then calling GMT modules directly using clib.Session.call_module() calls.
liked it better than calling high level pygmt.Figure() and its methods
this way I don’t need to learn pygmt (or GMT.lj) equivalents of gmt command line arguments and options, and still have all the power of python available thanks to pygmt virtual files.
saved a lot of effort rewriting a python script in bash just because gmt clip is missing among pygmt.Figure() methods.