Showing posts with label gmt. Show all posts
Showing posts with label gmt. Show all posts

Thursday, December 10, 2009

create a grd file from DEM

After creating a DEM file using SRTM.py script, now say we need to convert it to a GMT .grd file.

1. convert the binary file to a grd file:
xyz2grd [binary file] -G[output_grd_file] -Zh -I[pixel_size] -R[west]/[south]/[east]/[north]r
- don't forget the r at the end.
- in the -Zh, h is for 2 byte integer. change if needed
2. Sample to a different pixel size and sub-region:
grdsample [input_grd_file] -G[output_grd_file] -I[new_pixel_size] -R[west]/[south]/[east]/[north]r
- if a known number of pixels is desired, say in east-west, use: east = west+[needed_pixel_number - 1]*[new_pixel_size]
3. put zeros instead of null values -9999:
grdmath [input_grd_file] -9999 GT [input_grd_file] MUL = [output_grd_file]
- note that if you sample some of the null pixels will be with different values so use say -500 instead of -9999 as a threshold.

Friday, October 2, 2009

Reading GMT .GRD/NetCDF file using python

When reading GTM file in GRD format you actually read a NetCDF format.
using python and pylab to read it enable matplotlib to plot it.
say we have a grd file named z.grd which contains 3 variables: x,y and z.
reading z values in python is done by the following three lines:

from pylab import *
from scipy.io import netcdf_file as netcdf
data = netcdf('z.grd','r').variables['z'][::-1]

Now data is a numpy array containing z values in the file dimensions.
for more info look at: http://gfesuite.noaa.gov/developer/netCDFPythonInterface.html