Gridleme ve filtreleme

Sayfanın 15:13, 18 Şubat 2010 tarihindeki hali

• batitri verisini gridleme ve spike'lari maskeleme

 
 
### dosya ismi gir
set dosya = bathy.dat
 
gmtset MEASURE_UNIT inch
 
########get minmax coordinates
set R = `minmax $dosya -C | awk '{print int($1/10)*10-20,int($2/10)*10+20,int($3/10)*10-20,int($4/10)*10+20}'`
 
 
#Math:    aliased to set !:1 = `echo "!:3-$" | bc -l `
#Math x = $R[2]- $R[1]
#Math y = $R[4]- $R[3]
#echo $x $y
 
########### butun noktalari gridlemek icin
# gridding with 0.5 metre pixel size
xyz2grd -R$R[1]/$R[2]/$R[3]/$R[4] -I0.5 -F -V -Gtmp05.grd $dosya
# 1 metre with 1 metre pixel size
xyz2grd -R$R[1]/$R[2]/$R[3]/$R[4] -I1 -F -V -Gtmp1.grd $dosya
# gridding with 2 metre pixel size
xyz2grd -R$R[1]/$R[2]/$R[3]/$R[4] -I2 -F -V -Gtmp2.grd $dosya
# gridding with 5 metre pixel size
xyz2grd -R$R[1]/$R[2]/$R[3]/$R[4] -I5 -F -V -Gtmp5.grd $dosya
 
 
 
######### plot each to see the difference
#grd2cpt tmp05.grd -Crainbow -V -Z > ! palet.cpt
#grdimage tmp05.grd -Cpalet.cpt -R$R[1]/$R[2]/$R[3]/$R[4] -Jx0.007 -P -B200 > ! x.ps
#xv x.ps
 
#grd2cpt tmp1.grd -Crainbow -V -Z > ! palet.cpt
#grdimage tmp1.grd -Cpalet.cpt -R$R[1]/$R[2]/$R[3]/$R[4] -Jx0.007 -P -B200 > ! x.ps
#xv x.ps
 
#grd2cpt tmp2.grd -Crainbow -V -Z > ! palet.cpt
#grdimage tmp2.grd -Cpalet.cpt -R$R[1]/$R[2]/$R[3]/$R[4] -Jx0.007 -P -B200 > ! x.ps
#xv x.ps
 
#grd2cpt tmp5.grd -Crainbow -V -Z > ! palet.cpt
#grdimage tmp5.grd -Cpalet.cpt -R$R[1]/$R[2]/$R[3]/$R[4] -Jx0.007 -P -B200 > ! x.ps
#xv x.ps
 
 
#### resample and make all with the same pixel size
grdsample tmp1.grd -I0.5 -R$R[1]/$R[2]/$R[3]/$R[4] -Gtmp1S.grd
grdsample tmp2.grd -I0.5 -R$R[1]/$R[2]/$R[3]/$R[4] -Gtmp2S.grd
grdsample tmp5.grd -I0.5 -R$R[1]/$R[2]/$R[3]/$R[4] -Gtmp5S.grd
 
grd2cpt tmp5S.grd -Crainbow -V -Z > ! palet.cpt
grdimage tmp5S.grd -Cpalet.cpt -R$R[1]/$R[2]/$R[3]/$R[4] -Jx0.08 -P -B200 > ! x.ps
xv x.ps
 
##### stack them with highest resolution on top
## fill the gaps in the data with 0.5 m pixel size
## from the with 1 m, 2m and then 5 m 
grdmath tmp05.grd tmp1S.grd AND = tmpA.grd
grdmath tmpA.grd tmp2S.grd AND = tmpB.grd
grdmath tmpB.grd tmp5S.grd AND = spiky.grd
 
 
######### plot the stacked data with 0.5 m pixel size. 
### create an intensity (shade) file with sun angle of 315 (NW)
grdgradient spiky.grd -GspikyI.grd -Nt -A315
grd2cpt spiky.grd -Crainbow -V -Z > ! palet.cpt
grdimage spiky.grd  -Cpalet.cpt -R$R[1]/$R[2]/$R[3]/$R[4] -Jx0.07 -P -B200 > ! x.ps
xv x.ps &
 
### 3D plotting see the spikes !!!!!
echo $R
set nam = 3d.ps
set ranges = -R$R[1]/$R[2]/$R[3]/$R[4]  
set depth = ` grdinfo -C spiky.grd | awk '{print $6,$7}'`
grdview   spiky.grd -E40/20 $ranges -Jx0.08 -Cpalet.cpt -IspikyI.grd -Jz0.08 -Qi300 -V   > ! $nam
 
xv $nam
 
############# remove the spikes now
### make a shade and a slope file
grdgradient spiky.grd -Gintensity.grd -Dc -Sslope.grd -V
 
## convert pixels with a slope value higer than 1.2 to NaN (Not a Value)
grdclip slope.grd -Gmask.grd -Sa1.2/NaN -V
 
## mask spikes now
grdmath  spiky.grd mask.grd OR = nospike.grd
 
## have a look at it
grdgradient nospike.grd -GtmpI.grd -Nt -A315
grd2cpt nospike.grd -Crainbow -V -Z > ! palet.cpt
grdimage nospike.grd -ItmpI.grd -Cpalet.cpt -R$R[1]/$R[2]/$R[3]/$R[4] -Jx0.08 -P -B200 > ! x.ps
xv x.ps &
 
# make intensity file
grdgradient nospike.grd -GnospikeI.grd -Nt -A315
 
### 3D plotting 
set ranges = -R$R[1]/$R[2]/$R[3]/$R[4]  
 
grdview   nospike.grd -E10/20 $ranges -Jx0.08 -Cpalet.cpt -InospikeI.grd -Jz0.08 -Qi300 -V -Z0 -B250 > ! $dosya:r3d.ps
xv $dosya:r3d.ps &