Re: [ferret_users] SRTM Data with Ferret

[mike mahoney <m.mahoney@xxxxxxxxxxx>]

Ferreters, 

I forgot to CC this response to the list. It may be of
some use if others are interested in the SRTM data.


Hello Paulo

I have used the 3 arc-second SRTM data. I have attached a journal
file I use to convert the data to netcdf. 
I believe the format is the same as the SRTM30 which is a simple 16 bit
signed integer raster file of the data on a regular lat/long grid. 
The document "srtm30_documentation.pdf" (from the ftp site) gives
details.

There are some differences for which you will need to change the script.
The file name is made up from the lower left co-ordinates for the 3
arc-sec data (my script uses this) whereas the upper left corner is used
for the 30 arc-second files. The missing data value is -32768 for 3 a-s
and -9999 for 30 a-s and the file extension is "hgt" instead of "dem".

I have now modified the script (load_SRTM30.jnl) and tested it on two 
of the 30arc-second files and it seems to work.

Hope it is of some use

Cheers
mike

On Fri, 2008-01-11 at 16:47 -0300, Paulo Henrique wrote:
> Hi Ferreters,
> 
> Is there any way to use SRTM (Shuttle Radar Topography Mission) data
> with Ferret? 
> 
> More information at http://www2.jpl.nasa.gov/srtm/
> 
> The data is available at
> ftp://e0srp01u.ecs.nasa.gov/srtm/version2/SRTM30
> 
> Thanks in advance.
> 
> -- 
> Ms. Paulo Henrique Santiago de Maria
> Grupo de Modelagem Atmosférica 
> Departamento de Meteorologia e Oceanografia
> Fundação Cearense de Meteorologia e Recursos Hídricos
> Av. Rui Barbosa 1246 - CEP 60115-221
> Fortaleza, Ceará
> Fone: (85) 3101-1106 / 3101-1126


------------------------------------------------------------------------
 The information contained in this communication is  for the use of the 
 individual  or  entity  to  whom  it  is  addressed, and  may  contain 
 information which is the  subject of legal privilege and/or copyright. 
 If you have received this  communication in  error, please  notify the 
 sender by return E-Mail and delete the transmission, together with any 
 attachments, from your system. Thank you.
-------------------------------------------------------------------------

! Usage go load_SRTM $1 $2
!
! Arguments: $1 = integer data file to convert to netCDF
! Arguments: $2 = netcdf file to append to
! 	    


query/ignore $1"<First argument must be the integer HGT data file name"
query/ignore $2"<Second argument must be the file name for the netCDF file"
!query/ignore $3"<Third argument must be the title for the data set"

set data/save


!define axes and grid to read in upside down data

def axis/x=1:1201:1 xin ; def axis/y=1:1201:1 yin
def grid/x=xin/y=yin gin

! load the binary data
file/var=rawhgt/grid=gin/format=stream/type="i2"/swap "$1"

! use samplej to reverse the y grid
let yr=1202-y[g=gin]
let hgtr=samplej(rawhgt,yr)

! find the lower left corner co-ords from the filename
let LLLat=`substring("$1",2,2)`
let LLLon=`substring("$1",5,3)`


!define the geo axes and grid
let arc_secx3=(1E-03)/1.2
define axis/x=`LLLon`:`LLLon + 1.00`:`arc_secx3`/units=longitude longitude
define axis/y=`LLLat*-1`:`LLLat*-1 + 1.00`:`arc_secx3`/units=latitude latitude
define grid/Y=latitude/X=longitude geo_grid

!assign the rectified data to the geo grid
let hgt=hgtr[g=geo_grid,gx=@asn,gy=@asn]

set var/bad=-32768 hgt
set var/title="SRTM"/units="metres" hgt

save/file="$2"/append hgt

! Usage go load_SRTM30 $1 $2
!
! this version for the 30 arc-second SRTM data - untested **********
!
! Arguments: $1 = integer data file to convert to netCDF
! Arguments: $2 = netcdf file to append to
! 	    


query/ignore $1"<First argument must be the integer DEM data file name"
query/ignore $2"<Second argument must be the file name for the netCDF file"
!query/ignore $3"<Third argument must be the title for the data set"

set data/save


!define axes and grid to read in upside down data
! 30 arc-second grid is 50latx40lon degrees - 6000x4800 points

def axis/x=1:4800:1 xin ; def axis/y=1:6000:1 yin
def grid/x=xin/y=yin gin

! load the binary data
file/var=rawhgt/grid=gin/format=stream/type="i2"/swap "$1"

! use samplej to reverse the y grid
let yr=6001-y[g=gin]
let hgtr=samplej(rawhgt,yr)

! find the upper left corner co-ords from the filename
! This could come from the .hdr or .dmw files
let ULLat=`substring("$1",6,2)`
let ULLon=`substring("$1",2,3)`
def sym lat_tag=`substring(("$1"),5,1)`
def sym lon_tag=`substring(("$1"),1,1)`
if `"($lat_tag)" EQ "S"` then let ULLat=`ULLat * -1.0`
if `"($lon_tag)" EQ "W"` then let ULLon=`ULLon * -1.0`

!define the geo axes and grid
let arc_secx30=10/1200
define axis/x=`ULLon`:`ULLon + 40.00 - arc_secx30`:`arc_secx30`/units=longitude longitude
define axis/y=`ULLat`:`ULLat - 50.00 + arc_secx30`:`arc_secx30`/units=latitude latitude
define grid/Y=latitude/X=longitude geo_grid

!assign the rectified data to the geo grid
let hgt=hgtr[g=geo_grid,gx=@asn,gy=@asn]

set var/bad=-9999 hgt
set var/title="SRTM"/units="metres" hgt

save/file="$2"/append hgt