Re: BEDMAP with GMT
Sumant Jha <[email protected]> Tue, 24 Jun 2014 15:39:40 -0600
| Newsgroups | gmane.comp.gis.gmt.user |
|---|---|
| Message-ID | <[email protected]> |
Hi All,
Raising this topic again.
For some reason, the following solution does not works with the final release of bedmap2 which is available at https://secure.antarctica.ac.uk/data/bedmap2/
I have the old bm2_draft_bed.txt, and its output comes exactly as what Tom has attached. However, when using bedmap2_bed.txt (which I think, is the later version of bm2_draft_bed.txt), the output comes as attached.
Can anyone help me understand what has changed since the draft version to final version of BEDMAP2? Or is it that I am doing something wrong in my code? I am using same code that Tom has given in his email with minor changes like replacing the bm2_draft_bed.txt by bedmap2_bed.txt, using the bedmap2_bed.txt as input, and changing the -R to -180/180/-90/-60.
The code for bedmap2_bed and its output are attached.
#!/bin/sh
# bedmap_ant.sh
#
#
# Created by Sumant Jha on 12/27/12.
# Based on TOM JORDAN’s script, dated Nov 21, 2012 in [GMT-HELP] BEDMAP with GMT
# Modified on 06/24/2014
#
# script to get degree grid of BEDMAP2 data
# Note: For some reason, this script is only working on bm2_draft_bed.txt
# clean and order xyz data
echo "***********************************************"
echo “ Cleaning and Ordering X Y Z data in desired sequence"
awk ' NR>1 { print $1*1, $2*1, $3*1 } ' bedmap2_bed.txt | awk ' { print $2, $3, $1 } ' > bedmap
# projecting bedmap data to degrees (default bedmap comes in UTM) NOTE -C option
echo " "
echo " "
echo "***********************************************"
echo " Projecting BEDMAP2 to Degrees"
mapproject bedmap -I -R-180/180/-90/-60 -Js0/-90/-71/1:1 -C -F -V> temp1
# Make new 1 min grid
echo " "
echo " "
echo "***********************************************"
echo " Creating 1 min grid for BEDMAP2"
surface temp1 -R-180/180/-90/-60 -I1m -V -Gbedmap_1min.grd
# make colour pallet
echo " "
echo " "
echo "***********************************************"
echo " Creating color pallett for displaying data"
grd2cpt -Z -V bedmap_1min.grd > bedmap1.cpt
gmtset PAPER_MEDIA A4
# Draw Antarctica Bedmap grid
echo " "
echo " "
echo "***********************************************"
echo "Drawing Antarctica BEDMAP grid"
grdimage bedmap_1min.grd -Cbedmap1.cpt -R-180/180/-90/-60 -Js0/-90/-71/1:50000000 -B30g300/15g15wEsN -K -V -P> BEDMAP_ant.ps
# Draw coast
pscoast -W1 -R -Js -O >> BEDMAP_ant.ps
# Create contours spaced at 200 intervals (optional)
#grdcontour Bedmap_ant.grd -J -R -W1/80/80/80 -C200 -P -S4 -O -A1000+f3+k80/80/80+s8t -G2i/10 -Djunk -V -K >> BEDMAP_Ross.ps
echo " "
echo " "
echo "***********************************************"
echo "File Generated"
#View
gs BEDMAP_ant.ps
On Nov 21, 2012, at 6:46 AM, Jordan, Tom A.R.M. <[email protected]> wrote:
> Hi All
>
> One solution to the problem is below.
>
> First get the XYZ data for BEDMAP2 from the same ftp site Roi mentioned (it’s in a different directory).
> ftp://ftp.bas.ac.uk/ptf/bm2_ftp/text_version_of_bed/
>
> You can then re-project the data to degrees with mapproject, the trick is the –C option.
>
> Re grid the degree data and off you go.
>
> Alternatively you could dump the coast line from pscoast (-m option).
> Project the coast line to stereographic meters, same map project but no –I.
> Then simply plot the projected coast on top of the bedmap grid extracted as Roi suggested using psxy –Jx.
>
> Hope this is useful.
>
> Tom
>
> # script to get degree grid of BEDMAP2 data, its slow but works
> # clean and order xyz data
> awk ' NR>1 { print $1*1, $2*1, $3*1 } ' bm2_draft_bed.txt | awk ' { print $2, $3, $1 } ' > temp
>
> # project data to degrees NOTE -C option
> mapproject temp -I -R0/360/-90/0 -Js0/-90/-71/1:1 -C -F > temp1
>
> # Make new 5 min grid
> surface temp1 -R0/360/-90/-60 -I5m -Gbedmap_5min.grd
>
> # make colour pallet
> grd2cpt -Z bedmap_5min.grd > temp.cpt
>
> # draw grid
> grdimage bedmap_5min.grd -Ctemp.cpt -R0/360/-90/-60 -Js0/-90/-71/1:50000000 -K > temp.ps
>
> # draw coast
> pscoast -W1 -R -Js -O >> temp.ps
>
> # view
> gv temp.ps
>
>
> From: World-Wide GMT Usage and Help Mailing List [mailto:[email protected]] On Behalf Of Roi Granot
> Sent: 20 November 2012 17:49
> To: [email protected]
> Subject: Re: [GMT-HELP] BEDMAP with GMT
>
> There is a new Bedmap compilation (Bedmap2, doi:10.5194/tcd-6-4305-2012) which you can download here: ftp://ftp.nerc-bas.ac.uk/ptf/bm2_ftp/readme%20and%20data%20files/
>
> With gdal you can convert it:
> dal_translate -of GMT bm2_draft_bed.tif bedmap2.grd
>
> and plot:
> grdimage -Jx-.000001/-.000001 -Ibedmap2.grd.int -B5000000 -R-3333500/3333500/-3333500/3333500 bedmap2.grd -Cbed.cpt -V > bedmap2.ps
>
> Roi
>
> Dr. Roi Granot
> Department of Geological & Environmental Sciences
> Ben Gurion University of the Negev, Beer Sheva
> P.O.Box 653 Israel 84105
> Tel: 972-8-6477-509
>
>
>
>
>
> On Nov 20, 2012, at 3:43 PM, Marcus Badger wrote:
>
>
> Hi Sumant,
>
> I downloaded the "Grounded ice-sheet bed elevation" in Ascii format from BAS from here:
>
> http://www.antarctica.ac.uk//bas_research/data/access/bedmap/download/
>
> Then used xyz2grd to convert to grid. If you already have the grid then it might be fine. What do you get when you run
>
> grdinfo BEDMAP_groundbed.grd
>
> For my grid I get:
> groundbed.grd: Title: groundbed.grd
> groundbed.grd: Command: xyz2grd groundbed.asc -E -Ggroundbed.grd
> groundbed.grd: Remark:
> groundbed.grd: Pixel node registration used
> groundbed.grd: Grid file format: nf (# 18) GMT netCDF format (float) (COARDS-compliant) [DEFAULT]
> groundbed.grd: x_min: -2661600 x_max: 2678400 x_inc: 5000 name: x nx: 1068
> groundbed.grd: y_min: -2149967 y_max: 2195033 y_inc: 5000 name: y ny: 869
> groundbed.grd: z_min: -2495.36010742 z_max: 4364.85400391 name: z
> groundbed.grd: scale_factor: 1 add_offset: 0
>
> So you can see that the x and y are not in latlon but metres relative to the south pole so in my script I just set -R to the min and max. As this is already projected data you should then just need to overlay the appropriate latlon grid using psbasemap -Js with the appropriate parameter to match the projection of the grid. The docs give "we used the Polar Stereographic projection with 71°S as the latitude of true scale and 0°E as the central meridian" although I'm not sure what parameters you therefore need for -Js
>
> Marcus.
>
>
> ----------------------------------------------------------
> Dr Marcus Badger
> Lecturer
>
> School of Earth Sciences
> University of Bristol
> Wills Memorial Building
> Queen's Road
> BRISTOL BS8 1RJ
> U.K.
>
> Office: G40
> Tel: +44 (0) 117 95 45491 or +44 (0) 117 33 16795
> Mobile: +44 (0) 774 07 96811
> Email: [email protected]
> ----------------------------------------------------------
>
> On 16/11/2012 17:20, Sumant Jha wrote:
>
> Hi Markus,
>
> Pardon my ignorance on this, as I am very new to GMT. Do we need to convert the BEDMAP grid to some other format? Where are you getting the xyz data for the same?
>
> I was using the BEDMAP grid that came with Geoware.
>
> You are right about the -R values. I noticed the same but was not sure, how to go about converting my -R values to that of BEDMAP values. Your version is something what I am trying to achieve, and then overlay the lat/longs. Not sure how to do that.
>
> Sumant
>
>
> On Nov 15, 2012, at 2:41 AM, Marcus Badger <[email protected]> wrote:
>
> Hi Sumant,
>
> How far have you got? ie in what way is it not coming out right? I've managed to make this into a grid using xyz2grd and plot it using
>
> grdimage groundbed.grd -V -Cpalette.cpt -B2000000NseW -R-2661600/2678400/-2149967/2195033 -JX5i -Yc -Xc -P -K >Antarctica.ps
>
> but the grid is not in lat-lon but from the docs page "In this system, coordinates are expressed in metres with the origin at the South Pole." (which might be where you problems start) and so I am unsure how to convert this to something you can overlay coastline and gridlines. My version attached.
>
>
> Marcus.
>
> ----------------------------------------------------------
> Dr Marcus Badger
> Lecturer
>
> School of Earth Sciences
> University of Bristol
> Wills Memorial Building
> Queen's Road
> BRISTOL BS8 1RJ
> U.K.
>
> Office: G40
> Tel: +44 (0) 117 95 45491 or +44 (0) 117 33 16795
> Mobile: +44 (0) 774 07 96811
> Email: [email protected]
> ----------------------------------------------------------
>
> On 13/11/2012 20:34, Sumant Jha wrote:
> Guys - I fear this question has been asked a lot before - but, i cannot
> seem to find an answer to it. I am trying to use the BEDMAP data in GMT.
> The problem it seems to me is the projection. I tried using Jx and then
> Js to overlay the lat-longs - but, somehow, it did not come out right. I
> am hoping if someone out here have worked with GMT and BEDMAP data - and
> could may be help me solve this issue.
> I am using the following code to plot it:
>
> grdimage
> /Users/sumantjha/Documents/Geoware/DATA/misc/BEDMAP_groundbed.grd -V
> -C/Users/sumantjha/Documents/Work/GMT/GMT4.5.8/share/cpt/GMT_gebco.cpt
> -B20NseW -R155/180/-80/-75 -Jx10i -X1.2i -Y6i -K -V -P > Antarctica.ps
>
> pscoast -R -Js0/-90/50i/1 -W -Df -Ir/.5/114/187/243 -G30/144/255 -V -O
> -K -P >> Antarctica.ps
>
> Any help will be appreciated.
>
> Thanks
> Sumant
> To unsubscribe, send the message "signoff gmt-help" to
> [email protected]
>
> To unsubscribe, send the message "signoff gmt-help" to [email protected]
> <Antarctica.jpg>
>
> To unsubscribe, send the message "signoff gmt-help" to [email protected]
>
>
> To unsubscribe, send the message "signoff gmt-help" to [email protected]
>
>
> To unsubscribe, send the message "signoff gmt-help" to [email protected]
>
>
>
> --
> This message (and any attachments) is for the recipient only. NERC
> is subject to the Freedom of Information Act 2000 and the contents
> of this email and any reply you make may be disclosed by NERC unless
> it is exempt from release under the Act. Any material supplied to
> NERC may be stored in an electronic records management system.
> To unsubscribe, send the message "signoff gmt-help" to [email protected] <temp.jpg>
Mailing list for GMT discussions of all kinds. If you are not sure you have found a bug, discuss it here first.
To formally report bugs or request features, please register and add New Issue on gmt.soest.hawaii.edu
To unsubscribe, send the message "signoff gmt-help" to [email protected]
Note: gmt-help will become obsolete on Sept 1, 2014 - please use forum on gmt.soest.hawaii.edu instead.