Convert Mercator grid to WGS84 deg
Nico Fournier <[email protected]>
| Newsgroups | gmane.comp.gis.gmt.user |
|---|---|
| Message-ID | <[email protected]> |
Dear all, I have a New Zealand bathymetry grid from NIWA in Mercator and would like to both 1. plot it directly for quick perusal and 2. convert it to geographical coordinates (WGS84 deg) for later use. 1. is very mildly successful while 2. is a no goer so far. Could anyone help me with any of those or - in an ideal world - both? note: I am using gmt4.5.9 on Mac OS X 10.9. On the NIWA website, the info about the dataset is: #------------------------------ Projection Mercator (WGS84 Datum) at 41°S true latitude and 100°E central meridian (EPSG 3994). CANZ, 2008: New Zealand Region Bathymetry, 1:4 000 000, 2nd Edition. NIWA Chart, Miscellaneous Series No. 85. (http://www.niwa.co.nz/our-science/oceans/bathymetry/further-information) #------------------------------ Please note that I downsampled the original NIWA grid to 2.5km to reduce the size of the test grid. Running grdinfo gives the following: #------------------------------ RegionalBathymetry2008_2.5km_v2.grd: Pixel node registration used RegionalBathymetry2008_2.5km_v2.grd: Grid file format: nf (# 18) GMT netCDF format (float) (COARDS-compliant) [DEFAULT] RegionalBathymetry2008_2.5km_v2.grd: x_min: 4795705.5515 x_max: 7824455.5515 x_inc: 2498.96864686 name: x nx: 1212 RegionalBathymetry2008_2.5km_v2.grd: y_min: -5915595.45788 y_max: -2067845.45788 y_inc: 2500.16244314 name: y ny: 1539 RegionalBathymetry2008_2.5km_v2.grd: z_min: -10045.1953125 z_max: 3338.73242188 name: z RegionalBathymetry2008_2.5km_v2.grd: scale_factor: 1 add_offset: 0 #------------------------------ 1. plotting: I can plot the map in cartesian system using the following. Is it the most appropriate way to do it though since it seems to distort the map? #------------------------------ ingrid='RegionalBathymetry2008_2.5km_v2.grd' area=`grdinfo -I- $ingrid` proj=-JX5i psbasemap $proj $area -B500000::WES -K -P >> $outpsfile grdgradient $ingrid -Gingrid.inten -A300/30 -Nt grd2cpt $ingrid -S-10000/4000/100/ -Z -Cocean > ingrid.cpt grdimage $ingrid -Cingrid.cpt $area $proj -Iingrid.inten -P -O -K >> $outpsfile psscale -D2.5i/5.4i/5i/0.2ih -S -Cingrid.cpt -B2000:"Elevation a.s.l. (m)":/:: -A -O -P >> $outpsfile #------------------------------ 2. conversion to geographical grid My understanding so far is that this can be achieved using grdproject. My feeble attempt is: >> grdproject $ingrid -G$outgrid -Jm100/-41/1:4000000 -I -R<xmin>/<xmax><ymin>/<ymax> Problem #1 is that I do not know the extent of the map in WGS84 deg for -R. Problem #2 is that I may have misunderstood the whole grdproject altogether (which is a distinct possibility…); in which case problem #1 may be irrelevant. Would anyone please have some pointers? Cheers Nico ------------------------------------------------- Nico Fournier, Senior Volcano Geophysicist GNS Science - Te Pu Ao, Wairakei Research Centre Taupo, New Zealand Web: www.gns.cri.nz - www.geonet.org.nz - www.volcanoloco.org Email: [email protected] ”Knowledge is knowing that a tomato is a fruit; wisdom is knowing not to put it in a fruit salad.”- Brian O’Driscoll Notice: This email and any attachments are confidential. If received in error please destroy and immediately notify us. Do not copy or disclose the contents. To unsubscribe, send the message "signoff gmt-help" to [email protected]