Plotting only in projected cartesian coordinates

Q&A
1

I have data in cartesian coordinates (easting northing from Ordnance Survey projection) and would like to plot only using this system. However, I’m struggling to understand what projection command to provide in order that I can show coastline data.
Here’s my code so far:

fig = pygmt.Figure()
fig.grdimage(
    "@earth_relief_01s",
    region=region_deg,
    projection=projection_gmt,
    shading="+a45+nt0.7",
    cmap="gray",
)
fig.coast(shorelines=True, region=region_deg)

maxabs = vd.maxabs(Bouguer_Anomaly)

pygmt.makecpt(cmap="polar", series=[-maxabs, maxabs])
fig.plot(x=long, y=lat, fill=Bouguer_Anomaly, cmap=True, style="c2p", projection=projection_gmt)
fig.basemap(projection=projection_gmt, frame=True)
fig.colorbar(cmap=True, frame=["a50f25", "x+lBouguer Anomaly", "y+lmGal"])
fig.plot(data=np.array([[ed,sd,wd,nd]]), style='r+s', pen="3p,yellow")
fig.show()

where ed,sd,wd,nd are the east, south, west and north values for an area, Bouguer_Anomaly is the data I’m looking to plot and region_deg is the region I’m plotting.

Note that I’m currently plotting using projection_gmt = "T-2/49/15c" since this is the projection of OSGB36 coordinates. However, I’d instead like to use easting northing data to plot, and therefore would like to define the conversion (i.e where the origin of the coordinates is, in order that pygmt knows what coast data to plot underneath my Bouguer_Anomaly data.

The main purpose behind this is that I want to reduce error accrued from constantly converting between cartesian and geographic coordinates (with the analysis processes I’m using requiring cartesian data).

I’m assuming the way forward would be to use pygmt.grdproject, but I’m still confused as to how to define the projection later on in the figure.

The system I’m using has true origin easting 40000, northing -10000; true lat long origin 49ºN,2ºW, scale factor on the central meridian of 0.9996012717, and uses the Transverse Mercator projection with the Airy 1830 ellipsoid.

Lastly, I should note that I’m using Verde, a gridding python library, meaning that vd.maxabs simply chooses the maximum value from the Bouguer_Anomaly data, for colorbar scaling purposes.

2

as your data are already projected, you can plot your data using Cartesian projection: projection = X6c or projection=x6c changing size or scale as needed.

To overlay GMT coastlines you need to correctly specify your projection to pscoast:

This may or may not work, see more below.

You would need to indicate your projection ellipsoid as shown here PROJ_ELLIPSOID
This is done using pygmt.config:

pygmt.config(PROJ_ELLIPSOID="6378137")

Check GMT manual for the Transverse Mercator projection projection, indicated by letter T or t: 6.3.2. Transverse Mercator projection (-Jt -JT). However I can see no way of specifying origin easting and northing, these true origin easting 40000, northing -10000 using (py)gmt -R syntax.

There is another possibility: to save coastlines in the region of interest from coast to files and reproject these to your coordinate system using mapproject. For reprojection you can use more or less any source and destination coordinate system using e.g. PROJ specifications, while the subset that actually works for plotting and drawing map frames is limited which is the problem for directly using pygmt.coast. The reprojected coastlines are free from this issue and can be potted using pygmt.plot() and the same Cartesian projection using projection=X6c.

Or, as you guessed, you could reproject your Cartesian grid to geographic coordinates and do all your plotting in geographical coordinates. Specifying the source ellipsoid and projection would also be necessary.