### 2015/04/28 Yuta Ando (Mie University)

#!/bin/csh

#setting for GMT
gmtset ANNOT_FONT_SIZE_PRIMARY = 10p
gmtset ANNOT_FONT_PRIMARY = Helvetica
gmtset DEGREE_SYMBOL = ring
gmtset PLOT_DEGREE_FORMAT = ddd:mm:ssF


#setting for map
set range=135:40/138:20/32:20/35:0 #136:20/137:20/34:20/35:10   #range of map
set range2=129/146/30/46     #Sub map
set size=18                   #size of map
set size2=10
set reso=f                    #resolution for map(c<l<i<h<f)
set xanot=a20mf10mg5m              #scale for horizontal axis
set yanot=a20mf10mg5m            #scale for vertical axis
set xanot2=a240mf120m
set yanot2=a240mf120m
set zanot=a2000f500g500
set view=150/30

set in0=ETOPO1_Ice_c_gmt4
set out=map                   #output file name


wget http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/cell_registered/netcdf/ETOPO1_Ice_c_gmt4.grd.gz
gunzip -d ETOPO1_Ice_c_gmt4.grd.gz
xyz2grd etopo1_ice_c_f4.flt -Getopo1.grd -I1m/1m -F -R-180/180/-90/90 -ZTLh -V

#make color palet
makecpt -V -Cglobe > globe2.cpt

#draw map
#color shading
grdimage ${in0}.grd -Cglobe2.cpt -JM${size} -R${range} -B${xanot}/${yanot}WSne -V -P -K > ${out}.ps
#grdview ${in0}.grd -Qs -Wc1 -JZ${size} -Cglobe2.cpt -E${view} -JM${size} -R${range}/-8000/8000 -B${xanot}/${yanot}/${zanot}:"Height [m]":WSneZ+ -V -P -K > ${out}.ps

#contours
grdcontour ${in0}.grd -JM -C500 -A500 -W2 -L-8000/-100 -R -B -V -P -K -O >> ${out}.ps
grdcontour ${in0}.grd -JM -C10 -A20 -W2 -L-50/-1 -R -B -V -P -K -O >> ${out}.ps
#grdcontour ${in0}.grd -E${view} -JZ -JM -C500 -A500 -W2 -L-8000/-1 -R -B -V -P -K -O >> ${out}.ps
grdcontour ${in0}.grd -JM -C500 -A500 -W2 -L1/8000 -R -B -V -P -K -O >> ${out}.ps
#grdcontour ${in0}.grd -E${view} -JZ -JM -C500 -A500 -W2 -L1/8000 -R -B -V -P -K -O >> ${out}.ps

#color scale
psscale -Cglobe2.cpt -D4c/-1.2c/7c/0.5ch -B10000/:m: -P -V -K -O >> ${out}.ps
psbasemap -R -JM -Lf136:10/32:33/20/20+l+jt -V -P -K -O >> ${out}.ps
#psbasemap -E${view} -R -JZ -JM -Lf138/32:30/33/50+l+jt -K -O >> ${out}.ps

#coast
pscoast -R -JM -D${reso} -W2 -P -V -K -O >> ${out}.ps
#pscoast -E${view} -R -JZ -JM -D${reso} -W3 -P -V -O >> ${out}.ps

#Sub map
grdimage ${in0}.grd -Cglobe2.cpt -JQ${size2} -R${range2} -B${xanot2}/${yanot2}wsne -X7.5 -Y0.5 -V -P -O >> ${out}.ps

#convert postscript to png
ps2raster ${out}.ps -A -Tg
rm ${out}.ps 

#display map with GGV
eog ${out}.png

