# additional lines to pygmt lecture

## import extra libraries

import rioxarray
import xarray as xr
import pygmt

## load velocity from previous (InSAR) lecture(s)

inpath = 'https://gws-access.jasmin.ac.uk/public/nceo_geohazards/shared/temp/earmla/022D_04826_121209.vel.geo.tif'
vel = rioxarray.open_rasterio(inpath)
vel = vel.load()

## a little tweak the loaded data

vel = vel.squeeze('band')
vel = vel.drop('band')

## create a basic map

fig = pygmt.Figure()
pygmt.makecpt(cmap="roma", series=[-30,60], reverse=True)
fig.grdimage(grid=vel, cmap=True)
fig.basemap(frame="a")
fig.colorbar(frame="a5")
fig.show()