SRTM çizimi

Sayfanın 22:44, 15 Haziran 2009 tarihindeki hali

 
########################################################################
## SRTM veri indirme, GMT grid'e dönüştürme, birleştirme ve plot etme ##
########################################################################
 
# SRTM verileri aşağıdaki ftp sitesinden indirilebilir
 
ftp 152.61.4.9 
# kullanıcı anonymous 
# şifre email adresi
 
# /srtm/Documentation klasörünün altında Continent_def.gif isimli dosyayı indirip ilgilenen bölge 
# hangi kategoriye giriyor öğren.
 
cd  /srtm/version1
get Continent_def.gif
 
# Farzedelim Eurasia kategorisi olsun. Bu durumda 90 m'lik yani 3 ark saniyelik SRTM verinizin 
# bulunduğu yer şurasıdır:
 
cd /srtm/version2/SRTM3/Eurasia 
 
# burada 1x1 derecelik çerçevelere bölünmüş durumda srtm verileri. Her bir çerçevenin nereyi 
# kapsadığı isminden anlaşılabilir. Çerçeve isimlerinde o çerçevenin sol alt köşe koordinatı 
# verilmektedir. Örneğin 40-41 kuzey enlemleri ve 29-30 boylamları arasında kalan bölgenin 
# verisini içeren dosyanın adı N40E029.hgt.zip'dir. Bu çerçeve Bursa-Gemlik civarını kapsamaktadır. 
 
# Şimdi bunu ve yandaki çerçeveyi indirelim.
 
get N40E029.hgt.zip
get N40E028.hgt.zip
 
# zipleri açalım
 
unzip N40E029.hgt.zip
unzip N40E028.hgt.zip
 

 
# srtm verilerini indirip birleştirmenin en kolay yolu getsrtm3 
# [http://www.anandarooproy.com/gmt/getsrtm.tgz] komutunu kullanmaktır. 
 
# Yukarıdaki iki dosya şu şekilde de indirilebilir.
 
getsrtm3 eurasia -R28/30/40/41
 
# İndirme işleminden sonra final.grd isimli bir dosya üretilir bu iki çerçevenin 
# birleştirilmiş gmt grididir. Ancak grdinfo ile mutlaka bakınız ve yüksekik 
# değerlerinin doğru olduğundan emin olunuz. Değil ise aşağıda nasıl düzeltileceğine 
# dair kısma bakınız
 
 

SRTM verileri 2 byte binary formatındadır. Yani her bir yükseklik verisi için diskte ayrılan yer 2 byte'dır. Bu durumda bir SRTM çerçevesinin ne kadarlık bir yer kapladığını hesap edebiliriz

1 enlem derece (degree) = 60 dakika (minute) = 3600 saniye (second)'dir

İndirdiğimiz SRTM 3 arc saniye (1 arc saniyelik DEM'ler sadece Kuzey Amerika için dağıtılıyor). Her bir piksel 3 saniyelik bir uzunluğa sahip olduğuna göre bir derecelik mesafe için gerekli olan piksel sayısı 3600/3 = 1200'dür.

Dolayısıyla, her bir çerçevede 1x1 derecelik bir alan olduğuna göre piksel sayımız 1200x1200 olması gerekiyor. Ancak piksel registration değil de grid registration yapıldığı için birer piksel fazla oluyor (yani 1201x1201).

Peki niye SRTM 90 m diyoruz???

Çünkü 1 derece = 111 km olduğuna göre 1 arc saniye = 111000/3600 = ~30m yapar. Buradan 3 saniye de 3x30m = 90 m (aslında tam olarak 92.5 m) yapar

Peki 90x90 kare şeklinde midir? Hayır!

Çünkü 3 saniye, eğer ekvator civarında değilseniz, kuzey-güney yönünde tam olarak aynı mesafeye takabül etmez . Enlemler arasındaki mesafe değişmez; dünyanın etrafını yaklaşık 40000 km kabul edersek 1 enlem derecenin km cinsinden değişimi 40000/360 ~= 111 km'dir. Boylamlar arasındaki mesafe ise bulunan enleme göre değişir. Ekvatordan kutuplara doğru gidildikçe mesafe azalır. Dünyanın yarıçapının 6367449 metre olduğu kabul edildiğinde X enlemi boyunca boylamlar arası mesafe için şu formül kullanılır:

 mesafe = 3.14/180 * cos(x) * 6367449 

Marmara bölgesinde enlem 40 derece olarak alınırsa bu mesafenin yaklaşık 86 km olduğu ortaya çıkar. Buradan da doğu batı yönünde piksel büyüklüğü 86/1200~=71 m olduğu çıkar. Dolayısıyla bu bölgede pikseller 71x92 m ebatındadır.

 GTOPO30 ve SRTM30 30 saniye aralıklıdır dolayısıyla 30x30 = ~900 m, yani SRTM90'nın 10 katı

 
# Her bir piksel 2 byte'lık bir yer tuttuğuna göre çerçevenin boyutunun 1201x1201x2=2884802 byte 
# olması gerekir. Bunu doğrulamak için: 
 
ls -la
 
# yapınca aşağıdaki gibi bir çıktı görürsünüz
 
# -rw-r--r--  1 ziyadin  staff  2884802 May 27  2005 N40E029.hgt
 
# ve buradan da dosya boyutunun gerçekten de  2884802 olduğu görülür. Aksi türlü bir yanlışlık var demektir.
 
# Şimdi bu dosyaları xyz2grd komutu ile gridleyip birleştirelim.
 
# -R ile DEM'in sınırları verilir: -R/lonMin/lonMax/latMin/latMax
# -I ile derece cinsinden piksel büyüklüğü: 1/1200 = 0.000833333
# -Zh ile verinin 2 byte binary olduğu belirtilir
# -G ile output grd dosya ismi verilir
 
xyz2grd N40E029.hgt -R29/30/40/41 -I0.00083333 -Zh  -GN40E029.grd
xyz2grd N40E028.hgt -R28/29/40/41 -I0.00083333 -Zh  -GN40E028.grd
 
# Şimdi grdinfo programı  ile yükseklik değerlerine bir bakalım mantıklı mı
 
grdinfo N40E029.grd
 
#  deyince şu satırlar gözükür
#   N40E029.grd: Title: N40E029.grd
#   N40E029.grd: Command: xyz2grd -R29/30/40/41 -I0.00083333 N40E029.hgt -Zh -GN40E029.grd
#   N40E029.grd: Remark: 
#   N40E029.grd: Gridline node registration used
#   N40E029.grd: Grid file format: nf (# 18)
#   N40E029.grd: x_min: 29 x_max: 30 x_inc: 0.000833333 name: x nx: 1201
#   N40E029.grd: y_min: 40 y_max: 41 y_inc: 0.000833333 name: y ny: 1201
#   N40E029.grd: z_min: -32768 z_max: 32521 name: z
#   N40E029.grd: scale_factor: 1 add_offset: 0
 
# Görüldüğü gibi bende z_min ve z_max değerleri -32768 ve 32521
# Burada bir hata var!!!
 
# Sizde de böyle çıkıyorsa o zaman bunun düzeltilmesi gerekiyor.
# Bunun sebebi binary dosyada byte-order'ın  (big veya little endian) yanlış olmasıdır. 
# Byte order'u swap etmek için -Z  opsiyonuna w eklenmesi gerekir. 
 
xyz2grd N40E029.hgt -R29/30/40/41 -I0.00083333 -Zhw  -GN40E029.grd
xyz2grd N40E028.hgt -R28/29/40/41 -I0.00083333 -Zhw  -GN40E028.grd
 
# Yeniden bakalım düzelmiş mi
 
grdinfo N40E029.grd
 
# dersek şunları görürüz
#   .....
#   N40E029.grd: z_min: -32768 z_max: 2531 name: z
#   ....
 
# x_max düzelmiş ama z_min halen kötü. Bu değeri sıfırlamamız gerekiyor. 
# Sıfırlamak için grdclip komutunu kullanabiliriz.
 
# -Sb0/NaN opsiyonu ile 0 den küçük tüm değerler [[NaN]] (Not A Number) olarak atanır ve 
# yeni grid -G opsiyonunda verilir. Deniz seviyesinin altına olan yerler var malum bu durumda 
# 0 yerine -1000 gibi bir değer verilebilir.
 
grdclip N40E029.grd -Sb0/NaN -GN40E029C.grd
grdclip N40E028.grd -Sb0/NaN -GN40E028C.grd
 
# Yeniden bakalım
 
grdinfo N40E029C.grd
 
#dersek
#
# ....
#N40E029C.grd: z_min: 0 z_max: 2531 name: z
# ....
 
grdinfo N40E028C.grd
#....
#N40E028C.grd: z_min: 0 z_max: 1251 name: z
#...
 
# bunlar yan yana olduğu için yine GMT grdpaste ile birleştirilip -G opsiyonu ile birleştirilen 
# dosyaya mozaiksrtm.grd ismi verilebilir.
 
#########################  birleştirme #################################
grdpaste N40E029C.grd N40E028C.grd -Gmozaiksrtm.grd
 
grdinfo mozaiksrtm.grd
 
#  mozaiksrtm.grd: Title: mozaiksrtm.grd
#  mozaiksrtm.grd: Command: grdpaste N40E029C.grd N40E028C.grd -Gmozaiksrtm.grd
#  mozaiksrtm.grd: Remark: 
#  mozaiksrtm.grd: Gridline node registration used
#  mozaiksrtm.grd: Grid file format: nf (# 18)
#  mozaiksrtm.grd: x_min: 28 x_max: 30 x_inc: 0.000833333 name: x nx: 2401
#  mozaiksrtm.grd: y_min: 40 y_max: 41 y_inc: 0.000833333 name: y ny: 1201
#  mozaiksrtm.grd: z_min: 0 z_max: 2531 name: z
#  mozaiksrtm.grd: scale_factor: 1 add_offset: 0
 
# Görüldüğü gibi doğu batı yönünde piksel adeti ikiye katlandı.
 
 
 
 
#################### Şimdi bunu plot edip görelim #####################
# bir grd değişkeni tanımlayıp tek tek her seferinde yazmaktan kurtulalım
set grd = mozaiksrtm.grd
 
# Renkli palet yapalım
grd2cpt -Crainbow -Z $grd > ! color.cpt
 
# Gölge dosyası oluşturalım
grdgradient  $grd -A315 -Nt -Gshade.grd
 
# Çizim için derece rakam formatı gibi bazı gmt ayarlarına ihtiyacımız var: 
gmtset PLOT_DEGREE_FORMAT +D MEASURE_UNIT inch
 
# alanı ve projeksiyonu (JM = mercator projeksiyon) R değişkenine atayalım:
set R = "-R28/30/40/41 -JM6"
 
# dosya ismi verelim
set name = marmararenkli.ps
 
psbasemap $R -B0.2 -K > ! $name
grdimage $R -Ccolor.cpt  $grd -Ishade.grd -O -K  >> $name
pscoast -Df -O $R -W -Ia >> $name 
xv $name &
 
 
# gri ton olsun ama  tonlama yükseltiye göre değişmesin; gölgeye göre değişşin.
# Palet oluşturalım:
# Kolonlar: altyükseklik R G B üstyükseklik R G B
# R,G ve B: 0-255 arası RGB skalası renk değerleri. Örneğin kırmızı: 255 0 0
echo -10000 180 180 180 10000 180 180 180 > ! gray.cpt
 
set name = marmaragri.ps
psbasemap $R -B0.2 -K -P> ! $name
grdimage $R -Cgray.cpt $grd -Ishade.grd -O -K>> $name
pscoast -Df -S120 -W -O -R -JM -Ia/2/blue -K >> $name
psbasemap -Lf29.8/40.1/29.8/40.1/20 -R -JM -O -B.2  >> $name
xv $name &
 

--Ziyadin 22:30, 9 Nisan 2009 (CEST)