SRTM GEBCO mozaik
From GMT Türkiye Wiki
######################################################################## ## 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/Documentation get Continent_def.gif # Farzedelim Eurasia katagorisi olsun. Bu durumda 90 m'lik yani 3 ark saniyelik SRTM verinizin # bulunduğu yer şudur cd /srt/version2/SRTM3/Eurasia # burada 1x1 derecelik çerçevelere bölünmüş durumda srtm verileri. Herbir çerçevenin nereyi # kapsadığı isminden anlaşılabilir. Çerçeve isminlerinde o çerçevenin sol alt köşe kordinatı # verilmektedir. Örneğin 40-41 kuzey enlemleri ve 29-30 boylamları arasında kalan bölgenin # veririni içeren dosyanın adı N40E029.hgt.zip'dir. Bu çeç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 verileri 2 byte binary formatındadır. Yani her bir yüksekik verisi için ayrılan yer 2 byte'dır. # Bu durumda her bir çerçevenin ne kadarlık bir yer kapladığını hesap edebiliriz. # 1 derece (degree) = 60 dakika (minute) = 3600 saniye (second)'dir # indirdiğimiz SRTM 3 arc saniye ( 1 arc saniye verisi sadece kuzey amerika için var sadece) # 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 # Dolayısıyla her bir çerçevede 1x1 derecelik bir alan olduğuna göre piksel sayımız 1200x1200 # olması gerekiyor. Ancak grid değilde [[piksel kaydı]] 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 = ~30 m yapar. Buradan # 3 saniye de 3x30 = 90 m yapar # GTOPO30 ve SRTM30 30 saniye aralıklıdır dolayısıyla 30x30 = ~900 m, yani srtm 90'nin 10 katı # Herbir piksel 2 byte lik bir yer tutuğuna göre çerçevenin boyutu 1201x1201x2=2884802 byte # olması gerekir. Bunu doğrulamak için ls -la # yapınca aşağıdaki gibi bir çıktı görürüsünüz # -rw-r--r-- 1 ziyadin staff 2884802 May 27 2005 N40E029.hgt # ve buradanda boyutun gerçektende 2884802 olduğu anlaşı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!!! # Sizdede böyle çıkıyorsa o zaman bunun düzeltilmesi gerekiyor. # Bunun sebebi binary dosyada byte-order'ın (big veya little endiean) yanlış olmasıdır. # 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 olduğ yerler var malum bu durumda # 0 yerine -1000 felan 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 yanyana 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 tektek 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 # bazı gmt ayarları gmtset PLOT_DEGREE_FORMAT +D MEASURE_UNIT inch # dosya ismi verelim set name = marmara.ps # alanı ve projeksiyonu (JM = markatör) setleyelim 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 #gri ton olsun ama yükseltiye göre değişmesin tonlama, gölgeye göre değişşin. # 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 & ############################################ ## GEBCO veri indirme srtm ile mozaikleme ## ############################################ # Web arayüyü kullanarak istenilen diktörtgen bir alan seçilip veya kordinatları verilip indirilebilir. # Default olarak indirilen dosya aslında hazır GMT grd'si ancak uzantısı nc olarak veriliyor her # nedense. Bu nedenle kafa karışmaması için soyadını grd'ye çevirmeniz yeterli. # Tüm grid de indirilebilir vey 1 GB'ın üstünde. Bu dosya grdinfo yaparsak şunu görebiliriz # (bendeki dosya halen nc soyadı taşıyor :( grdinfo GEBCO_08.nc #yaparsak şunları görüyoruz # GEBCO_08.nc: Title: GEBCO_08 Grid # GEBCO_08.nc: Command: 20090202 # GEBCO_08.nc: Remark: # GEBCO_08.nc: Pixel node registration used # GEBCO_08.nc: Grid file format: cs (# 8) # GEBCO_08.nc: x_min: -180 x_max: 180 x_inc: 0.00833333 name: user_x_unit nx: 43200 # GEBCO_08.nc: y_min: -90 y_max: 90 y_inc: 0.00833333 name: user_y_unit ny: 21600 # GEBCO_08.nc: z_min: -10977 z_max: 8685 name: user_z_unit # GEBCO_08.nc: scale_factor: 1 add_offset: 0 # Görüldüğü üzer bu grid tüm dünyayı kapsıyor. # Tüm grid indirilmişşe kesilmesi gerekiyor ve bunun için de grdcut komutu kullanılabilir. # 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 # Gebco verisinde piksel aralığı 30 saniye yani SRTM in 3 katı. Dolayısıyla bunu SRTM in # piksel boyutuna düşürelim grdsample komtu ile. grdsample -I0.000833333 Gebco_cut.grd -Ggebco_cutS.grd -T # grdinfo gebco_cutS.grd yaparsak aşağıda görüldüğü gibi piksel büyüklüğü (0.000833333), # alanı (28/30/40/41) ve dolayısıyla piksel sayısı (nx ve ny) aynı olur SRTM mozaiksrtm.grd ile. # gebco_cutS.grd: Title: gebco_cutS.grd # gebco_cutS.grd: Command: grdsample -I0.000833333 Gebco_cut.grd -Ggebco_cutS.grd -T # gebco_cutS.grd: Remark: # gebco_cutS.grd: Gridline node registration used # gebco_cutS.grd: Grid file format: nf (# 18) # gebco_cutS.grd: x_min: 28 x_max: 30 x_inc: 0.000833333 name: longitude [degrees_east] nx: 2401 # gebco_cutS.grd: y_min: 40 y_max: 41 y_inc: 0.000833333 name: latitude [degrees_north] ny: 1201 # gebco_cutS.grd: z_min: -1250.12 z_max: 2442.18 name: z # gebco_cutS.grd: scale_factor: 1 add_offset: 0 # 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 # Şimdi iki veri seti mozaiklenebilir çünkü tıpa tıp aynı iki grd. Bunun için Gebco gridinin kara # alanlarını NaN yapıp grdmath komutunu kullanabiliriz. # 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 #bu mask uygulayalım gebco verisine grdmath mask.grd gebco_cutS.grd ADD = gebmasked.grd # srtm ile maskelenmiş gebco verisini mozaikleyelim grdmath gebmasked.grd mozaiksrtm.grd AND = mozaik.grd #Geçmiş olsun mozaik hazırdır !!!!! #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
--Ziyadin 19:53, 1 Nisan 2009 (CEST)