#! /bin/csh -f
#
# Special purpose hackery to convert a "img grid" to geographic coordinates.
#
# NOTE: the right most column and top row will be lost, pad input but this much.
#
# The output of img2xyz (and img2web) is geographical units, but there are "extra rows". 
# These extra rows are  result of well known mercator projection distorting in y direction.
#
# What to do?
#
# interpolate x/y from geographic/mercator to geo/geo at twice the resolution in x and y
# this will effectively regrid the mercator to geo/geo. The reason we don't surface
# at the original x_inc, is that the mercator y nodes are not same as in geo y. 
# But they are almost the same for latitudes below about 80. 
# 
# That should have been trivial but we'll have to hack this grid to get what we want
# which is a PIXEL grid. The reason for that is "surface", is a NODE in / NODE out
# utility. 
# 
# IOW it gives NODE out for PIXEL input. Which is wrong.
# We could just toggle a normal grid with grd edit, 
# except in this case our grid is mercator. 
#
# So we'll have to play some games!

set srcFile = 'lame.xyz'
set inputRegion = '-R/0/240/-60/180'
set inputResolutionSec = 60

set outputResolutionSec = `echo $inputResolutionSec | awk '{ printf "%.16f\n", $1/2 }'`
echo $outputResolutionSec
#FIXME these should be something like -I30c ; but for debugging we skip the "c"
set inputResolution = `echo $inputResolutionSec | awk '{ printf "-I%s\n", $1 }'`
set outputResolution = `echo $inputResolutionSec | awk '{ printf "-I%.16f\n", $1/2 }'`

set grdSampleOptions = '-Ql'

surface   -V $srcFile             -G$srcFile.surface.grd $inputRegion $outputResolution
grdsample -V $srcFile.surface.grd -G$srcFile.surface.grdsample.grd -F $outputResolution $grdSampleOptions 

# At this point we have converted MERCATOR PIXEL xyz to GEO PIXEL grd...

# debug
minmax $srcFile $srcFile.surface.grd2xyz.xyz
grdinfo   -V $srcFile.surface.grd $srcFile.surface.grdsample.grd -L2 
grd2xyz   -V $srcFile.surface.grd           >! $srcFile.surface.grd2xyz.xyz
grd2xyz   -V $srcFile.surface.grdsample.grd >! $srcFile.surface.grdsample.grd2xyz.xyz

sort -n -k2 -k1 $srcFile
sort -n -k2 -k1 $srcFile.surface.grd2xyz.xyz
sort -n -k2 -k1 $srcFile.surface.grdcut.grd2xyz.xyz
sort -n -k2 -k1 $srcFile.surface.grdsample.grd2xyz.xyz

# Flipping surface output from NODE to PIXEL will change -R .
# Account for this: shift grid left and down by half of x_inc/y_inc
# 
# set hackedRegion = '-R/-15/165/-75/105'
# grdedit -V $srcFile.surface.grd $hackedRegion
# grdinfo -V -L2 $srcFile.surface.grd
# 
# Flipping surface output from NODE to PIXEL will change -R, but we just hacked that.
# In effect we have surfaced PIXEL in to PIXEL out; except we now have extrapolated stuff
# on the left & bottom edges, and lost one row & column on top/right 
# 
# grdedit -V -T  $srcFile.surface.grd 
# grdinfo -V -L2 $srcFile.surface.grd
# 
# Delete extrapolated stuff on left & bottom. The data on the top & right is lost.
# This is why we needed a pad on the input region.
# 
# FIXME
# grdcut  -V $srcFile.surface.grd -G$srcFile.surface.grdcut.grd $inputRegion
# grdinfo -V -L2 $srcFile.surface.grdcut.grd
# grd2xyz -V $srcFile.surface.grdcut.grd >! $srcFile.surface.grdcut.grd2xyz.xyz

exit




# set inputRegion = '-R/0/240/-60/180'
# set blockMedianFiltered = $srcFile.blockmedian.xyz
# 
# minmax $srcFile
# blockmedian -V $srcFile $inputRegion $outputResolution >! $blockMedianFiltered
# minmax $blockMedianFiltered
# sort -n -k2 -k1 lame.xyz.blockmedian.xyz > ! foo.xyz
# diff -w foo.xyz lame.xyz
# 
# 
# surface -V $srcFile -G$blockMedianFiltered.surface.grd $inputRegion $outputResolution
# grdinfo -V -L2 $blockMedianFiltered.surface.grd
# grd2xyz -V $blockMedianFiltered.surface.grd >! $blockMedianFiltered.surface.grd2xyz.xyz
# 
# grdsample -V  $blockMedianFiltered.surface.grd -G$blockMedianFiltered.surface.grdsample.grd -F $outputResolution -Ql
# grdinfo -V -L2 $blockMedianFiltered.surface.grdsample.grd
# grd2xyz -V $blockMedianFiltered.surface.grdsample.grd >! $blockMedianFiltered.surface.grdsample.grd2xyz.xyz
# 
# grdedit -V $blockMedianFiltered.surface.grd $hackedRegion
# grdinfo -V -L2 $blockMedianFiltered.surface.grd
# 
# # Flipping surface output from NODE to PIXEL will change -R, but we just hacked that.
# # In effect we have surfaced PIXEL in to PIXEL out; except we now have extrapolated stuff
# # on the left & bottom edges, and lost one row & column on top/right 
# 
# grdedit -V -T  $blockMedianFiltered.surface.grd 
# grdinfo -V -L2 $blockMedianFiltered.surface.grd
# 
# # Delete extrapolated stuff on left & bottom. The data on the top & right is lost.
# # This is why we needed a pad on the input region.
# 
# grdcut  -V $blockMedianFiltered.surface.grd -G$blockMedianFiltered.surface.grdcut.grd $inputRegion
# # grdinfo -V -L2 $blockMedianFiltered.surface.grdcut.grd
# 
# At this point we have converted MERCATOR PIXEL xyz to GEO PIXEL grd...
# debug
# grdinfo -L2 $blockMedianFiltered*.grd
# 
# grd2xyz -V $blockMedianFiltered.surface.grdcut.grd >! $blockMedianFiltered.surface.grdcut.grd2xyz.xyz
# minmax $blockMedianFiltered
# minmax $blockMedianFiltered*xyz
# sort -n -k2 -k1 $blockMedianFiltered
# sort -n -k2 -k1 $blockMedianFiltered.surface.grd2xyz.xyz
# sort -n -k2 -k1 $blockMedianFiltered.surface.grd2xyz.xyz
# sort -n -k2 -k1 $blockMedianFiltered.surface.grdcut.grd2xyz.xyz
#