Gridding PACES gravity data with GMT Stuart Wier Nov 13 2007 UNAVCO wier @ unavco .org 1. Quick Start 2. More details 1. Quick start: install GMT; see http://gmt.soest.hawaii.edu/; note that this installs NetCDF too. move columns from original data file, columns nos 5,3,and 12 (long, lat, and bouger anomaly value) to cols 1,2,3 in new "xyz" file cat paces_wyoming_gravity.txt | perl -ne 'chomp;@s=split(/ +/,$_);print "@s[5],@s[3],@s[12]\n";' > paces_wyoming_gravity_bouger.xyz remove 'bad' top lines in paces_wyoming_gravity_bouger.xyz run this GMT command: blockmean -R-112/-104/41/45 -I30c -V paces_wyoming_gravity_bouger.xyz > paces_wyoming_gravity_bouger_blockmean.xyz and this command: surface paces_wyoming_gravity_bouger_blockmean.xyz -R-112/-104/41/45 -I30c -Gpaces_wyoming_gravity_Bouger.grd -V This makes a NetCDF file of a sort, but without enough metadata for the IDV to understand it. To use this NetCDf file in the IDV you need to have identify variables (x,y) as latitude and longitude. Convert the binary netcdf file to its ascii CDL equivalent with the NetCDF program command: ncdump paces_colorado_gravity_Bouger_values.grd > paces_colorado_gravity_Bouger_values.cdl Edit the .cdl file and replace all instances of the variable name "x" with "lon" and all instances of the varaible name "y" with lat. You may also rename the field variable "z" to a more infromative name such as in this case "bouger_anomaly." and enter the "units= " lines as below. the top of the cdl file should look like: netcdf paces_wyoming_gravity_Bouger { dimensions: lon = 961 ; lat = 481 ; variables: float lon(lon) ; lon:long_name = "longitude" ; lon:units = "degrees_east" ; lon:actual_range = -112., -104. ; float lat(lat) ; lat:long_name = "latitude" ; lat:units = "degrees_north" ; lat:actual_range = 41., 45. ; float Bouger_anomaly(lat, lon) ; Bouger_anomaly:long_name = "Bouger_anomaly" ; Bouger_anomaly:actual_range = -293.153076171875, -83.6245193481445 ; Bouger_anomaly:coordinates = "lon lat" ; Now convert the new CDL file into a ".nc" NetCDF file usable in the IDV: with the netcdf command: ncgen -o paces_colorado_gravity_Bouger_values.nc paces_colorado_gravity_Bouger_values.cdl 2. More details get GeoNET gravity idata and put into files for IDV use get grav and mag data from http://paces.geo.utep.edu/welcome.shtml which has long-lat-data column order suitable for GMT esp from http://paces.geo.utep.edu/gdrp/Search.aspx for gravity or from http://gis.utep.edu/Projects/DataCatalog/tabid/57/Default.aspx The new terrain corrected United States gravity database is available through the GeoNet link below. The terrain corrections were calculated by Mike Webring of the U. S. Geological Survey using a digital elevation model and a technique based on the approach of Donald Plouff. The reduction of these data will be updated this year with modern geodetic datums and a higher precision digital elevation model. In the present version, latitude and longitude values are referenced to NAD27 (North American Datum 27; horizontal datum) and NAD83, elevation values in meter are referenced to NGVD29 (Vertical datum) and NGVD88. Both of these datums are listed in the dataset download. and thence esp. from http://irpsrvgis00.utep.edu/repositorywebsite/ "GeoNet Gravity and Magnetic Dataset Repository Pan American Center for Earth and Environmental Sciences (PACES)" cut out desired columns of data: for the old paces .dat data files get - the 9th column in this case, plus cols 1 and 2: cat paces_colorado_gravity.dat | perl -ne 'chomp;$col=9;$s=$col-3; /^([\w\-\.]+)\t([\w\-\.]+)\t([\w\-\. ]+\t){$s}([\w\-\.]+)/;print "$1,$2,$4\n";' > paces_colorado_gravity_Bouger_values.xzy for the new paces *.txt data files use cat paces_wyoming_gravity.txt | perl -ne 'chomp;@s=split(/ +/,$_);print "@s[3],@s[5],@s[12]\n";' > paces_wyoming_gravity_bouger.xyz to get lat (col 3) long (col 5), and bouger anom (col 12). the result *.xzy file is an ascii csv file; (long and lat in that order must be first output columns) remove the first 3 to 6 bad lines which are relicts from the source file header lines. run blockmean like: blockmean -R-110/-101/36/42 -I78c -V paces_colorado_gravity_Bouger_values.xzy > paces_colorado_gravity_Bouger_blockmean.xzy or blockmean -R-118/-103/40/49 -I30c -V paces_idaho_montana_wyoming_gravity_bouger.xyz > paces_idaho_montana_wyoming_gravity_bouger_blockmean.xyz or blockmean -R-112/-104/41/45 -I30c -V paces_wyoming_gravity_bouger.xyz > paces_wyoming_gravity_bouger_blockmean.xyz (you may have to delete a few lines in the source data file that have no data value) results like: blockmean: Given domain implies x_inc = 0.0216867 blockmean: Given domain implies y_inc = 0.0216606 blockmean: W: -110 E: -101 S: 36 N: 42 nx: 416 ny: 278 blockmean: Using a total of 3.53 Mb for all arrays. blockmean: Working on file paces_colorado_gravity_Bouger_values.xzy blockmean: Calculating block means blockmean: N read: 61775 N used: 61775 N outside_area: 0 N cells filled: 21613 grid it: surface paces_colorado_gravity_Bouger_blockmean.xzy -R-109/-102/37/41 -I30c -Gpaces_colorado_gravity_Bouger_values.grd -V note the -R should be the same area or smaller than in the blockmean command; the -I argument should be the same number or larger than in blockmean command. ??? set tension how? what value? results: surface: Given domain implies x_inc = 0.0216718 surface: Given domain implies y_inc = 0.0216216 surface: W: -109 E: -102 S: 37 N: 41 nx: 323 ny: 185 surface: WARNING: Your grid dimensions are mutually prime. surface: HINT: Choosing nx = 324, ny = 192 might cut run time by a factor of 107.9177 surface: HINT: Choosing nx = 360, ny = 192 might cut run time by a factor of 103.50099 surface: HINT: Choosing nx = 324, ny = 216 might cut run time by a factor of 101.5696 surface: HINT: Choosing nx = 360, ny = 200 might cut run time by a factor of 98.96357 surface: HINT: Choosing nx = 384, ny = 192 might cut run time by a factor of 98.139326 surface: HINT: Choosing nx = 360, ny = 216 might cut run time by a factor of 92.68345 surface: HINT: Choosing nx = 400, ny = 192 might cut run time by a factor of 91.257045 surface: HINT: Choosing nx = 400, ny = 200 might cut run time by a factor of 89.046667 surface: HINT: Choosing nx = 432, ny = 192 might cut run time by a factor of 87.091478 surface: HINT: Choosing nx = 324, ny = 240 might cut run time by a factor of 86.334161 surface: Minimum value of your dataset x,y,z at: surface: -106.361 39.2009 -339.29 surface: Maximum value of your dataset x,y,z at: surface: -102.517 40.9699 -107.975 surface: 496 unusable points were supplied; these will be ignored. surface: You should have pre-processed the data with block-mean, -median, or -mode. surface: Grid Mode Iteration Max Change Conv Limit Total Iterations surface: 1 D 250 0.169843 0.0357331 250 surface: Fit info: N data points N nodes mean error rms error curvature surface: 16416 60264 2.75099e-05 0.0202926 3230.9 The ".grd" file made is a NetCDF file. To use this NetCDf file in the IDV you need to have icentify which variable are latitude and longitude. Convert the binary netcdf file to its ascii CDL equivalent with ncdump paces_colorado_gravity_Bouger_values.grd > paces_colorado_gravity_Bouger_values.cdl Edit the .cdl file and replace all instances of the variable name "x" with "lon" and all instances of the varaible name "y" with lat. You may also rename the field variable "z" to a more infromative name such as in this case "bouger_anomaly." Now convert the new CDL file into a ".nc" NetCDF file usable in the IDV: ncgen -o paces_colorado_gravity_Bouger_values.nc paces_colorado_gravity_Bouger_values.cdl another way is to convert the output binary file paces_colorado_gravity_Bouger_values.grd to ascii column file: grd2xyz paces_colorado_gravity_Bouger_values.grd >! paces_colorado_gravity_Bouger_values_grid.xyz the result paces_colorado_gravity_Bouger_values_grid.xyz is space separated ascii file of columns the lat/long are in 3 sig fig; in htis example delta lat and long are 0.0216 degree there are 60264 lines, same a N nodes in surface output