#! /bin/csh -f

# img2xyz returns PIXEL registered, and clearly on a mercator grid.
#
# -AND- img2xyz return one extra column on the right and one extra row on the bottom.
#
# If you're fussy you'll have to either 
#	trim the xyz, 
#		ignore the xyz2grd error,
#			or hack the input
#
#img2web ../img/topo_14.1.img -11 41 179 221 1.0 >! w180n40.img2web.xyz

img2web ../img/topo_14.1.img -11 41 179 220.99165 1.0 > ! w180n40.img2web.xyz

minmax w180n40.img2web.xyz


#
# img2xyz outputs PIXEL grid, meaning the xy reported are  location of center of pixel
# img2xyz_80 outputs a 1m grid, we want 1m, but it's in mercator and we want geographic
#
# median filter onto a geographic pixel grid sounds good, but can't be correct
#

set img2xyzOutputName = 'w180n40.img2web.xyz'
set img2xyzOutputName = 'img2xyz'

set surfaceInputName  = $img2xyzOutputName
set surfaceOutputName = $img2xyzOutputName.surface

set grdeditOutputName = $surfaceOutputName
set grdeditInputName  = $grdeditOutputName.orig

# cp img2xyz_80.xyz $img2xyzOutputName.xyz
# sort -n -k3 $img2xyzOutputName.xyz >! $img2xyzOutputName.sorted.xyz

# Use surface to generate a NODE register geographical grid. We don't want NODE, but we 
# have no choice. The trick here is that to get the desired effect, we need to make sure 
# that we surface to a grid with a fine enough pitch that no two mercator nodes will be 
# in the same geographical node.
#
# We don't want to average, or median filter, we want to surface on a fine grid.
#
# Note this step is done at 30c to avoid two 1m mercator pixels going in one output pixel

set nodeRegisteredR = `jjSurfaceRegion_wrapper 179 221 -11 41 60 60`
surface $surfaceInputName.xyz -G$surfaceOutputName.grd $nodeRegisteredR -I30c= -V -VI
# FIXME: save a copy of what surface gave us for debugging
cp $grdeditOutputName.grd $grdeditInputName.grd
 
# convert to PIXEL registered
# NOTE: you can -NOT- combine the following two grdedit as that messes up WESN and inc...

grdedit $grdeditOutputName.grd -T -V
grdinfo -L2 $grdeditInputName.grd $grdeditOutputName.grd 

# convert the 0-360 lon to +/- 180
# NOTE that we need to use 30c on this step!!!

set nodeRegisteredR = `jjSurfaceRegion_wrapper -181 -139 -11 41 30 30`
grdedit $grdeditOutputName.grd $nodeRegisteredR -V
grdinfo -L2 $grdeditOutputName.grd 


exit

minmax $surfaceInputName.xyz 

# set blockmedianInputName = $surfaceOutputName
# set blockmedianOutputName = $blockmedianInputName.median
# 
# set xyz2grdInputName  = $blockmedianOutputName
# #set xyz2grdInputName  = $img2xyzOutputName
# set xyz2grdOutputName = $xyz2grdInputName




img2mercgrd -V -R/179/221/-11/41 ../img/topo_14.1.img -Gimg2grd.merc.grd -T1 -m1 -D
grdproject -V -Jm1 -R179/221/-11.0148512673/41.0050579259 -I img2grd.merc.grd -Gimg2grd.grd



minmax $surfaceInputName.xyz $blockmedianInputName.xyz $blockmedianOutputName.xyz
grdinfo -L0 $surfaceOutputName.grd $surfaceOutputName.edited.grd $surfaceOutputName.edited.pixel.grd $xyz2grdOutputName.grd $grdeditOutputName.grd


# grd2xyz $surfaceOutputName.edited.grd -R/-181/-139/-11/41 -S -V >! $blockmedianInputName.xyz
# blockmedian $blockmedianInputName.xyz -R/-181/-139/-11/41 -Q -I30c= $forcePixelReg -V >! $blockmedianOutputName.xyz
# sort -n -k3 $blockmedianOutputName.xyz >! $blockmedianOutputName.sorted.xyz
# 
# xyz2grd $xyz2grdInputName.xyz -G$xyz2grdOutputName.grd -I30c= -R/-181/-139/-11/41 $forcePixelReg -V
# 
# cp $grdeditInputName.grd $grdeditOutputName.grd 
# grdedit $grdeditOutputName.grd -R/-181/-139/-11/41 -V
# 

# now toggle grid from node to pixel registered
#
# this inevitably requires interpolation. But which is best -Ql, -Qn, or bicubic default?
# 
# grdsample $surfaceOutputName.edited.grd -T -G$surfaceOutputName.edited.pixel.grd
# grdsample $surfaceOutputName.edited.grd -T -Ganswer.grd -I1m
# 



echo 	img2grd $img -G$tmpFile1 -T1 -D $img2grdInc -V $region
cp /Volumes/publicSonarData/SRTM15/com/w180n40.img2xyz.grd $tmpFile2
	#remove mercator
# 	grdsample $tmpFile1 $region -I1m -G$tmpFile2 -Qn
	#resample to final resolution
grdsample $tmpFile2 -R/$w/$e/$s/$n $inc -G$dir/$stem.grd -Ql