SRTM GEBCO mozaik

From GMT Türkiye Wiki
(Sürümler arası farklar)
Jump to: navigation, search
 
1. satır: 1. satır:
 +
<csh>
 
#########################
 
#########################
 
## SRTM veri indirme, GMT grid'e dönüştürme ve plot etme
 
## SRTM veri indirme, GMT grid'e dönüştürme ve plot etme
266. satır: 267. satır:
 
pscoast -Df -O $R -W -Ia >> $name  
 
pscoast -Df -O $R -W -Ia >> $name  
 
xv $name
 
xv $name
 +
</csh>

Sayfanın 19:23, 1 Nisan 2009 tarihindeki hali

<csh>

    2. SRTM veri indirme, GMT grid'e dönüştürme ve plot etme
  1. SRTM verileri aşağıdaki ftp sitesinden indirilebilir

ftp 152.61.4.9

  1. kullanıcı anonymous
  2. şifre email adresi
  1. /srtm/Documentation klasörünün altında Continent_def.gif isimli dosyayı indirip ilgilenen bölge
  2. hangi kategoriye giriyor öğren.

cd /srtm/Documentation get Continent_def.gif

  1. Farzedelim Eurasia katagorisi olsun. Bu durumda 90 m'lik yani 3 ark saniyelik SRTM verinizin
  2. bulunduğu yer şudur

cd /srt/version2/SRTM3/Eurasia

  1. burada 1x1 derecelik çerçevelere bölünmüş durumda srtm verileri. Herbir çerçevenin nereyi
  2. kapsadığı isminden anlaşılabilir. Çerçeve isminlerinde o çerçevenin sol alt köşe kordinatı
  3. verilmektedir. Örneğin 40-41 kuzey enlemleri ve 29-30 boylamları arasında kalan bölgenin
  4. veririni içeren dosyanın adı N40E029.hgt.zip'dir. Bu çeçeve Bursa-Gemlik civarını kapsamaktadır.
  1. Şimdi bunu ve yandaki çerçeveyi indirelim.

get N40E029.hgt.zip get N40E028.hgt.zip

  1. zipleri açalım

unzip N40E029.hgt.zip unzip N40E028.hgt.zip


  1. SRTM verileri 2 byte binary formatındadır. Yani her bir yüksekik verisi için ayrılan yer 2 byte'dır.
  2. Bu durumda her bir çerçevenin ne kadarlık bir yer kapladığını hesap edebiliriz.
  1. 1 derece = 3600 saniye'dir
  1. indirdiğimiz SRTM 3 arc saniye ( 1 arc saniye verisi sadece kuzey amerika için var sadece)
  1. Her bir piksel 3 saniyelik bir uzunluğa sahip olduğuna göre bir derecelik mesafe için gerekli
  2. olan piksel sayısı 3600/3 = 1200
  3. Dolayısıyla her bir çerçevede 1x1 derecelik bir alan olduğuna göre piksel sayımız 1200x1200
  4. olması gerekiyor. Ancak grid değilde piksel kaydı yapıldığ için birer piksel fazla oluyor
  5. yani 1201x1201.
  1. Herbir piksel 2 byte lik bir yer tutuğuna göre çerçevenin boyutu 1201x1201x2=2884802 byte
  2. olması gerekir. Bunu doğrulamak için

ls -la

  1. yapınca aşağıdaki gibi bir çıktı görürüsünüz
  1. -rw-r--r-- 1 ziyadin staff 2884802 May 27 2005 N40E029.hgt
  1. ve buradanda boyutun gerçektende 2884802 olduğu anlaşılır. Aksi türlü bir yanlışlık var demektir.
  1. Şimdi bu dosyaları xyz2grd komutu ile gridleyip birleştirelim.
  1. -R ile dem'in sınırları verilir: -R/lonMin/lonMax/latMin/latMax
  2. -I ile derece cinsinden piksel büyüklüğü: 1/1200 = 0.000833333
  3. -Zh ile verinin 2 byte binary olduğu belirtilir
  4. -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

  1. Şimdi grdinfo programı ile yükseklik değerlerine bir bakalım mantıklımı

grdinfo N40E029.grd

  1. deyince şu satırlar gözükür
  2. N40E029.grd: Title: N40E029.grd
  3. N40E029.grd: Command: xyz2grd -R29/30/40/41 -I0.00083333 N40E029.hgt -Zh -GN40E029.grd
  4. N40E029.grd: Remark:
  5. N40E029.grd: Gridline node registration used
  6. N40E029.grd: Grid file format: nf (# 18)
  7. N40E029.grd: x_min: 29 x_max: 30 x_inc: 0.000833333 name: x nx: 1201
  8. N40E029.grd: y_min: 40 y_max: 41 y_inc: 0.000833333 name: y ny: 1201
  9. N40E029.grd: z_min: -32768 z_max: 32521 name: z
  10. N40E029.grd: scale_factor: 1 add_offset: 0
  1. Görüldğü gibi bende z_min ve z_max değerleri -32768 ve 32521
  2. Burada bir hata var!!!
  1. Sizdede böyle çıkıyorsa o zaman bunun düzeltilmesi gerekiyor.
  2. Bunun sebebi binary dosyada byte-order'ın (big veya little endiean) yanlış olmasıdır.
  3. 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

  1. Yeniden bakalım düzelmişmi

grdinfo N40E029.grd

  1. dersek şunları görürüz
  3. N40E029.grd: z_min: -32768 z_max: 2531 name: z
  1. x_max düzelmiş ama z_min halen kötü. Bu değeri sıfırlamamız gerekiyor.
  2. Sıfırlamak için grdclip komutunu kullanabiliriz.
  1. -Sb0/NaN opsiyonu ile 0 den küçük tüm değerler NaN (Not A Number) olarak atanır ve
  2. yeni grid -G opsiyonunda verilir. Deniz seviyesinin olduğ yerler var malum bu durumda
  3. 0 yerine -1000 felan verilebilir.

grdclip N40E029.grd -Sb0/NaN -GN40E029C.grd grdclip N40E028.grd -Sb0/NaN -GN40E028C.grd

  1. Yeniden bakalım

grdinfo N40E029C.grd

  1. dersek
  3. N40E029C.grd: z_min: 0 z_max: 2531 name: z

grdinfo N40E028C.grd

  1. .
  2. N40E028C.grd: z_min: 0 z_max: 1251 name: z
  1. bunlar yanyana olduğu için yine GMT grdpaste ile birleştirilip -G opsiyonu ile birleştirilen
  2. dosyaya mozaiksrtm.grd ismi verilebilir.

grdpaste N40E029C.grd N40E028C.grd -Gmozaiksrtm.grd

grdinfo mozaiksrtm.grd

  1. mozaiksrtm.grd: Title: mozaiksrtm.grd
  2. mozaiksrtm.grd: Command: grdpaste N40E029C.grd N40E028C.grd -Gmozaiksrtm.grd
  3. mozaiksrtm.grd: Remark:
  4. mozaiksrtm.grd: Gridline node registration used
  5. mozaiksrtm.grd: Grid file format: nf (# 18)
  6. mozaiksrtm.grd: x_min: 28 x_max: 30 x_inc: 0.000833333 name: x nx: 2401
  7. mozaiksrtm.grd: y_min: 40 y_max: 41 y_inc: 0.000833333 name: y ny: 1201
  8. mozaiksrtm.grd: z_min: 0 z_max: 2531 name: z
  9. mozaiksrtm.grd: scale_factor: 1 add_offset: 0
  1. Görüldüğü gibi doğu batı yönünde piksel adeti ikiye katlandı.
  1. Şimdi bunu plot edip görelim
  1. bir grd değişkeni tanımlayıp tektek her seferinde yazmaktan kurtulalım

set grd = mozaiksrtm.grd

  1. Renkli palet yapalım

grd2cpt -Crainbow -Z $grd > ! color.cpt

  1. Gölge dosyası oluşturalım

grdgradient $grd -A315 -Nt -Gshade.grd

  1. bazı gmt ayarları

gmtset PLOT_DEGREE_FORMAT +D MEASURE_UNIT inch

set name = marmara.ps set R = "-R28/30/40/41 -JM6" 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


  1. gri ton olsun ama yükseltiye göre değişmesin tonlama, gölgeye göre değişşin.
  2. Palet oluşturalım

echo -10000 180 180 180 10000 180 180 180 > ! gray.cpt

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 &


    2. GEBCO veri indirme srtm ile mozaikleme
  1. Web arayüyü kullanarak istenilen diktörtgen bir alan seçilip veya kordinatları verilip indirilebilir.
  2. Default olarak indirilen dosya aslında hazır GMT grd'si ancak uzantısı nc olarak veriliyor her
  3. nedense. Bu nedenle kafa karışmaması için soyadını grd'ye çevirmeniz yeterli.
  1. Tüm grid de indirilebilir vey 1 GB'ın üstünde. Bu dosya grdinfo yaparsak şunu görebiliriz

grdinfo GEBCO_08.nc

  1. GEBCO_08.nc: Title: GEBCO_08 Grid
  2. GEBCO_08.nc: Command: 20090202
  3. GEBCO_08.nc: Remark:
  4. GEBCO_08.nc: Pixel node registration used
  5. GEBCO_08.nc: Grid file format: cs (# 8)
  6. GEBCO_08.nc: x_min: -180 x_max: 180 x_inc: 0.00833333 name: user_x_unit nx: 43200
  7. GEBCO_08.nc: y_min: -90 y_max: 90 y_inc: 0.00833333 name: user_y_unit ny: 21600
  8. GEBCO_08.nc: z_min: -10977 z_max: 8685 name: user_z_unit
  9. GEBCO_08.nc: scale_factor: 1 add_offset: 0
  1. Görüldüğü gibi tüm dünya.


  1. Tüm grid indirilmişşe kesilmesi gerekiyor ve bunun için de grdcut komutu kullanılabilir.
  2. Gebco'yu yukarıda mozaiklediğimiz SRTM gridinin boyutunda keselim aşağıdaki komutla

grdcut -R28/30/40/41 GEBCO_08.nc -GGebco_cut.grd -fg

  1. Gebco verisinde piksel aralığı 30 saniye yani SRTM in 3 katı. Dolayısıyla bunu SRTM in
  2. piksel boyutuna düşürelim grdsample komtu ile.


grdsample -I0.000833333 Gebco_cut.grd -Ggebco_cutS.grd -T

  1. grdinfo gebco_cutS.grd yaparsak aşağıda görüldüğü gibi piksel büyüklüğü (0.000833333),
  2. alanı (28/30/40/41) ve dolayısıyla piksel sayısı (nx ve ny) aynı olur SRTM mozaiksrtm.grd ile.
  1. gebco_cutS.grd: Title: gebco_cutS.grd
  2. gebco_cutS.grd: Command: grdsample -I0.000833333 Gebco_cut.grd -Ggebco_cutS.grd -T
  3. gebco_cutS.grd: Remark:
  4. gebco_cutS.grd: Gridline node registration used
  5. gebco_cutS.grd: Grid file format: nf (# 18)
  6. gebco_cutS.grd: x_min: 28 x_max: 30 x_inc: 0.000833333 name: longitude [degrees_east] nx: 2401
  7. gebco_cutS.grd: y_min: 40 y_max: 41 y_inc: 0.000833333 name: latitude [degrees_north] ny: 1201
  8. gebco_cutS.grd: z_min: -1250.12 z_max: 2442.18 name: z
  9. gebco_cutS.grd: scale_factor: 1 add_offset: 0
  1. mozaiksrtm.grd: Title: mozaiksrtm.grd
  2. mozaiksrtm.grd: Command: grdpaste N40E029C.grd N40E028C.grd -Gmozaiksrtm.grd
  3. mozaiksrtm.grd: Remark:
  4. mozaiksrtm.grd: Gridline node registration used
  5. mozaiksrtm.grd: Grid file format: nf (# 18)
  6. mozaiksrtm.grd: x_min: 28 x_max: 30 x_inc: 0.000833333 name: x nx: 2401
  7. mozaiksrtm.grd: y_min: 40 y_max: 41 y_inc: 0.000833333 name: y ny: 1201
  8. mozaiksrtm.grd: z_min: 0 z_max: 2531 name: z
  9. mozaiksrtm.grd: scale_factor: 1 add_offset: 0
  1. Şimdi iki veri seti mozaiklenebilir çünkü tıpa tıp aynı iki grid. Bunun için Gebco gridinin kara
  2. alanlarını NaN yapıp grdmath komutunu kullanabiliriz.
  1. karaların NaN denizlerin 0 olduğu bir mask dosyası yapalım -N0/NaN ile

grdlandmask -Df gebco_cutS.grd -N0/NaN -Gmask.grd -R28/30/40/41 -I0.000833333

  1. bu mask uygulayalım gebco verisine

grdmath mask.grd gebco_cutS.grd ADD = gebmasked.grd

  1. srtm ile maskelenmiş gebco verisini mozaikleyelim

grdmath gebmasked.grd mozaiksrtm.grd AND = mozaik.grd

  1. Geçmiş olsun mozaik hazırdır !!!!!


  1. Plot et

gmtset PLOT_DEGREE_FORMAT +D MEASURE_UNIT inch set grd = mozaik.grd grd2cpt -Crainbow -Z $grd > ! color.cpt grdgradient $grd -A315 -Nt -Gshade.grd set name = marmozaik.ps set R = "-R28/30/40/41 -JM6" 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 </csh>

Retrieved from "https://atag.itu.edu.tr/w/index.php?title=SRTM_GEBCO_mozaik&oldid=1548"
Personal tools