Re: BEDMAP with GMT
Sumant Jha <[email protected]> Wed, 25 Jun 2014 09:02:13 -0600
| Newsgroups | gmane.comp.gis.gmt.user |
|---|---|
| Message-ID | <[email protected]> |
Thank you Paul and Tom. Tom: Thanks for clearing up my confusion about the two versions of BEDMAP2 data. I tried your xyz2grd command, but it gives me error: xyz2grd: Cannot open file –N-9999 Any idea what might be causing it? Paul: I am also trying your suggestions, and will report back on how it went. Regards Sumant On Jun 25, 2014, at 8:17 AM, Jordan, Tom A.R.M. <[email protected]> wrote: > Hi > > The original script works with bm2_draft_bed.txt as it is an ascii file with x, y, and z values. > > However, bedmap2_bed.txt is an ascii grid file. > > Try the following to convert the ascii grid file to a GMT format grid: > xyz2grd bedmap2_bed.txt -H6 -F -ZTLA –N-9999 -R-3333500/3333500/-3333500/3333500 -I1000 -GBEDMAP2_GMT.grd > > This should make you a GMT grid. You can then use grd2xyz to dump out an x, y and z ascii file just like bm2_draft_bed.txt, and run it through the original script, or follow Pauls suggestion below to avoid gridding over the pole. > > Tom > > From: Paul Wessel [mailto:[email protected]] > Sent: 25 June 2014 09:25 > To: [email protected] > Subject: Re: [GMT-HELP] BEDMAP with GMT > > I am surprised this ever worked. > Using surface to grid geographic data over a pole is a horrid idea since surface is best for Cartesian data or latitude-banded data. > If you absolutely need to convert this polar stereographic grid to degrees then I would do something like this: > 1. Create a blank grid with desired resolution for -R-180/180/-90/-60 > 2. Dump via grd2xyz and project coordinates with the -Js projection > 3. Sample the stereotraphic bedmap grid at the projected coordinates > 4. paste the result with the grd2xyz output and send via xyz2grd to make the desired grid > > Or something similar. > -p > > > On Jun 24, 2014, at 11:39 AM, Sumant Jha <[email protected]> wrote: > > > 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 > > <BEDMAP_ant.png> > > 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. > > 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. > > 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. > 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. 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.