# (29) Gridding spherical surface data using splines¶

Finally, we demonstrate how gridding on a spherical surface can be
accomplished using Green’s functions of surface splines, with or without
tension. Global gridding does not work particularly well in Cartesian
coordinates hence the chosen approach. We use
greenspline to produce a crude
topography grid for Mars based on radii estimates from the Mariner 9 and
Viking Orbiter spacecrafts. This data comes from *Smith and Zuber*
[Science, 1996] and is used here as a small (*N* = 370) data set we can
use to demonstrate spherical surface gridding. Since
greenspline must solve a *N* by
*N* matrix system your system memory may impose limits on how large data
sets you can handle; also note that the spherical surface spline in
tension is particularly slow to compute.

Our script must first estimate the ellipsoidal shape of Mars from the
parameters given by *Smith and Zuber* so that we can remove this
reference surface from the gridded radii. We run the gridding twice:
First with no tension using *Parker*’s [1990] method and then with
tension using the *Wessel and Becker* [2008] method. The grids are then
imaged with grdimage and
grdcontour and a color scale is placed between them.

```
#!/bin/bash
# GMT EXAMPLE 29
#
# Purpose: Illustrates spherical surface gridding with Green's function of splines
# GMT modules: makecpt, grdcontour, grdimage, grdmath greenspline, psscale, pstext
# Unix progs: rm, echo
#
ps=example_29.ps
# This example uses 370 radio occultation data for Mars to grid the topography.
# Data and information from Smith, D. E., and M. T. Zuber (1996), The shape of
# Mars and the topographic signature of the hemispheric dichotomy, Science, 271, 184-187.
# Make Mars PROJ_ELLIPSOID given their three best-fitting axes:
a=3399.472
b=3394.329
c=3376.502
gmt grdmath -Rg -I4 -r X COSD $a DIV DUP MUL X SIND $b DIV DUP MUL ADD Y COSD DUP MUL MUL Y \
SIND $c DIV DUP MUL ADD SQRT INV = PROJ_ELLIPSOID.nc
# Do both Parker and Wessel/Becker solutions (tension = 0.9975)
gmt greenspline -RPROJ_ELLIPSOID.nc mars370.txt -D4 -Sp -Gmars.nc
gmt greenspline -RPROJ_ELLIPSOID.nc mars370.txt -D4 -Sq0.9975 -Gmars2.nc
# Scale to km and remove PROJ_ELLIPSOID
gmt grdmath mars.nc 1000 DIV PROJ_ELLIPSOID.nc SUB = mars.nc
gmt grdmath mars2.nc 1000 DIV PROJ_ELLIPSOID.nc SUB = mars2.nc
gmt makecpt -Crainbow -T-7/15 > mars.cpt
gmt grdimage mars2.nc -I+ne0.75+a45 -Cmars.cpt -B30g30 -BWsne -JH0/7i -P -K -E200 \
--FONT_ANNOT_PRIMARY=12p -X0.75i > $ps
gmt grdcontour mars2.nc -J -O -K -C1 -A5 -Glz+/z- >> $ps
gmt psxy -Rg -J -O -K -Sc0.045i -Gblack mars370.txt >> $ps
echo "0 90 b)" | gmt pstext -R -J -O -K -N -D-3.5i/-0.2i -F+f14p,Helvetica-Bold+jLB >> $ps
gmt grdimage mars.nc -I+ne0.75+a45 -Cmars.cpt -B30g30 -BWsne -J -O -K -Y4.2i -E200 \
--FONT_ANNOT_PRIMARY=12p >> $ps
gmt grdcontour mars.nc -J -O -K -C1 -A5 -Glz+/z- >> $ps
gmt psxy -Rg -J -O -K -Sc0.045i -Gblack mars370.txt >> $ps
gmt psscale -Cmars.cpt -O -K -R -J -DJBC+o0/0.15i+w6i/0.1i -I --FONT_ANNOT_PRIMARY=12p -Bx2f1 -By+lkm >> $ps
echo "0 90 a)" | gmt pstext -R -J -O -N -D-3.5i/-0.2i -F+f14p,Helvetica-Bold+jLB >> $ps
# Clean up
rm -f *.nc mars.cpt
```