# Example 16 -- Grided bathymetry data, converted to image, then countoured # # Contributed 2001-May-14 by Peter S. Galbraith <GalbraithP@dfo-mpo.gc.ca> `Set image colorscale gebco' Set up colors for gebco charts of ocean/land altimetry. { open "awk 'BEGIN { \ for (i = 0; i < 256; i++) { \ if (i <= 278) {R=0.820; G=0.604; B=0.361;}\ if (i <= 232) {R=0.882; G=0.988; B=0.969;}\ if (i <= 227) {R=0.749; G=0.949; B=0.925;}\ if (i <= 223) {R=0.627; G=0.910; B=0.894;}\ if (i <= 209) {R=0.514; G=0.871; B=0.871;}\ if (i <= 185) {R=0.408; G=0.804; B=0.831;}\ if (i <= 139) {R=0.310; G=0.733; B=0.788;}\ if (i <= 93) {R=0.220; G=0.655; B=0.749;}\ if (i <= 46) {R=0.133; G=0.573; B=0.710;}\ if (i <= 0) {R=0.059; G=0.486; B=0.671;}\ print(R, G, B); \ } \ }' |" read image colorscale rgb close } set trace on set page landscape # Set my own variable to emphasize the Gulf only? Or the Ocean? .gulf_emphasisGSL. = 1 # Map lat/lon corners \xl = "73" \xr = "50" \yb = "43" \yt = "53" # Setup the map parameters set x format "%.0lf$\circ$W" set y format "%.0lf$\circ$N" set x type map W set y type map N set x name "" set y name "" set y margin 3 set x margin 3 set tic size {rpn ..tic_size.. 2 /} # Make smaller tic sizes (by half) set y size 14 if {rpn "/usr/local/bin/map-merc-scale-factor" file_exists} # Use my proj-based script to calc Mercator resizing if available \.factor. = system map-merc-scale-factor -a \yb -A \yt -o \xr -O \xl set x size {rpn ..ysize.. \.factor. *} else # else use Gri (Less than 0.2% difference anyway) resize x for maps end if set x axis \xl \xr -5 -1 # I want bigger tic marks at degree intervals, even if unlabelled set font size 0 set y axis \yb \yt 1 {rpn 1 6 /} draw y axis # Reset to label every second one to draw the labelled axes. set y axis \yb \yt 2 {rpn 1 6 /} set font to Helvetica set font size 12 draw axes # Get the x and y grid point locations from the pre-gridded `x y z' data file. if {rpn "y.grid" file_exists !} system zcat topex-gsl.xyz.gz | awk '{print $2}' | uniq > y.grid system zcat topex-gsl.xyz.gz | awk '{print $1}' | sort | uniq > x.grid end if open x.grid read grid x close open y.grid read grid y close # Generate the grid matrix if not already done if {rpn "topex-gsl.grid.gz" file_exists !} open topex-gsl.xyz.gz read columns x y z close set missing value 99999 # `x y z' data file is regularly-spaced, pre-gridded, # so use neighbour method convert columns to grid neighbor system rm -f topex-gsl.grid write grid to "topex-gsl.grid" system gzip --best topex-gsl.grid else open topex-gsl.grid.gz read grid data end if show grid if .gulf_emphasisGSL. Set image colorscale gebco set image range -500 50 set x type linear set y type linear set x format "%.0f m" draw image palette left -600 right 100 increment 100 \ box ..xmargin.. {rpn ..ymargin.. ..ysize.. + 2 + } \ {rpn ..xmargin.. ..xsize.. +} {rpn ..ymargin.. ..ysize.. + 3 + } else Set image colorscale gebco set image range -5000 500 set x type linear set y type linear set x format "%.0f m" draw image palette left -6000 right 1000 increment 1000 \ box ..xmargin.. {rpn ..ymargin.. ..ysize.. + 2 + } \ {rpn ..xmargin.. ..xsize.. +} {rpn ..ymargin.. ..ysize.. + 3 + } end if # Convert our `x y z' data grid to an image following above specs convert grid to image size 691 451 set clip postscript on # Clip in case we go over the edge by a wee bit draw image # Draw selected depth contours on the data grid using a thin line set line width .05 draw contour -50 unlabelled draw contour -100 -500 -100 unlabelled draw contour -1000 -3000 -1000 unlabelled # Draw some fillable coastline data set missing value -99 open gulfh_lewerIslands.lon_lat.gz read columns x y close read colornames from RGB "/usr/lib/X11/rgb.txt" set colour rgb 0.820 0.604 0.361 draw curve filled # Draw filled-in continents in this colour. set colour black draw curve # Draw again in black, without filling-in. draw axes frame open "/data/atlas/gmt/gulff-rivers-all.lon_lat" read columns x y set colour rgb {rpn 204 255 /} {rpn 221 255 /} {rpn 225 255 /} close draw curve open gulfh-level2.lon_lat.gz read columns x y close draw curve filled # Overlay lakes open gulfh-level3.lon_lat read columns x y close set colour rgb 0.820 0.604 0.361 draw curve filled # Overlay islands in lakes. set clip postscript off set colour black draw axes frame # Draw axes frame again to clean up quit