#! /bin/bash

PATH=${PATH}:/home/evan/bin
export PATH

# before running, make sure to tell the program where the netcdf file that contains your base topography is. In this example,
# I use GEBCO

topo_location="/home/evan/experiments/uplift_rate2/GEBCO/gebco_08_1m.grd"


# note that the grid spacing must be an integer! (in metres)

resolution=5000

min_ss=1000 # basal shear stress for areas outside of the domain


# create input files for the ice sheet program. Requires that you provide a CSV file of longitude/latitude values that is the boundary of your ice sheet that you want to create


outline_file_csv=$1 # long/lat file



echo 
echo "Entered convert_grid.sh using the following file:"
echo ${outline_file_csv}




# provide path to the long/lat shear stress grid file
shear_stress_ll=domains_ss.txt # created when create_shear_stress_file.sh is run



x_resolution=${resolution}
y_resolution=${resolution}

# Longitude region of interest, make sure it encompasses the provided outline at a minimum
minmax=$(awk -F',' '{if (NR == 1) {long_min = $1; long_max=$1; lat_min=$2; lat_max=$2} else { if ($1 > long_max) long_max=$1; if ($1 < long_min) long_min=$1; if ($2 > lat_max) lat_max=$2; if ($2 < lat_min) lat_min=$2;}} END{if (int(lat_max)+5 >= 90 ) lat_max=85; print int(long_min) - 5, int(long_max) + 5, int(lat_min)-5, int(lat_max)+5}' ${outline_file_csv})


echo "minmax, should be long (min max), lat (min max)"

echo ${minmax}

long_min=$( echo ${minmax} | awk '{print $1}')
long_max=$( echo ${minmax} | awk '{print $2}')

# latitude region of interest

lat_min=$( echo ${minmax} | awk '{print $3}')
lat_max=$( echo ${minmax} | awk '{print $4}')

echo ${long_min} ${long_max} ${lat_min} ${lat_max}

# Alberts equal area conic projection to convert from long/lat to Cartesian. 
# Obviously it is set up for NAICE, so the center is located near the center of the Keewatin dome at LGM

center_longitude=-108
center_latitude=61

southern_standard_parallel=47
northern_standard_parallel=74



map_width=20c # likely arbitrary


# this converts the outline file to be in cartesian coordinates

outline_file=outline.xyz

echo
echo "creating outline file"
echo


mapproject ${outline_file_csv}  -R${long_min}/${long_max}/${lat_min}/${lat_max} -F -C -Jb${center_longitude}/${center_latitude}/${southern_standard_parallel}/${northern_standard_parallel}/${map_width}  > ${outline_file}


xmin=$(sort -nk 1 ${outline_file} | head -n 1 | awk -v res=${resolution} '{print int(($1 -20*res) / res)*res}' )
echo xmin: ${xmin}
xmax=$(sort -nk 1 ${outline_file} | tail -n 1 | awk -v res=${resolution} '{print int(($1 +20*res) / res)*res}')
echo xmax: ${xmax}
ymin=$(sort -nk 2 ${outline_file} | head -n 1 | awk -v res=${resolution} '{print int(($2 -20*res) / res)*res}')
echo ymin: ${ymin}
ymax=$(sort -nk 2 ${outline_file} | tail -n 1 | awk -v res=${resolution} '{print int(($2 +20*res) / res)*res}')

echo "dimensions for ${outline_file}"
echo $xmin $xmax $ymin $ymax
echo
# first, make the basal shear stress file



execute_ss=y

if [ "${execute_ss}" = "y" ]
then

echo
echo "creating shear stress file"
echo

shear_stress_xyz_file=shear_stress.xyz

awk -v long_min=${long_min} -v long_max=${long_max} -v lat_min=${lat_min} -v lat_max=${lat_max} '{if ($1 >= long_min && $1 <= long_max && $2 >= lat_min && $2 <= lat_max) print $0 }' ${shear_stress_ll} > shear_stress_ll_temp

mapproject shear_stress_ll_temp  -R${long_min}/${long_max}/${lat_min}/${lat_max} -F -C -Jb${center_longitude}/${center_latitude}/${southern_standard_parallel}/${northern_standard_parallel}/${map_width}  > ${shear_stress_xyz_file}






blockmean -Re${xmin}/${xmax}/${ymin}/${ymax} -I${resolution} ${shear_stress_xyz_file} > bm.out

ss_grid="ss.grd"



triangulate bm.out -G${ss_grid} -I${resolution} -Re${xmin}/${xmax}/${ymin}/${ymax}



makecpt -Cwysiwyg -T0/120000/5000  -I > ss_shades.cpt
 grdimage ${ss_grid} -JX6i+ -P  -Css_shades.cpt > quick2.ps


num_x=$( ncdump -h ${ss_grid} | grep $'\tx = ' | sed -e 's/ //g' | sed 's/.*=//' | sed 's/;//g' )
num_y=$( ncdump -h ${ss_grid} | grep $'\ty = ' | sed -e 's/ //g' | sed 's/.*=//' | sed 's/;//g' )
echo dimensions of x and y
echo $num_x $num_y

ncdump  -v x  ${ss_grid}  | awk '/data:/{y=1;next}y' | sed -e 's/[x={}; ]//g' | tr -d '\n' | sed -e 's/,/\n/g' > x_values.txt

ncdump  -v y  ${ss_grid}  | awk '/data:/{y=1;next}y' | sed -e 's/[y={}; ]//g' | tr -d '\n' | sed -e 's/,/\n/g' > y_values.txt

xyz_file="index_ss.xyz"

# the -Cf option outputs it into array indices for reading into a Fortran program

grd2xyz ${ss_grid} -V -sa -Cf -N${min_ss} > ${xyz_file}

even_grid_file="even_grid_ss.dat"

nearest_int ${num_x} ${num_y} ${resolution} ${xyz_file} ${even_grid_file}

reduced_ss_file="reduced_ss.xyz"


reduce_dem ${even_grid_file} ${outline_file} ${reduced_ss_file}


binary_ss_file="ss.bin"

sed -e "s/${reduced_ss_file}/${binary_ss_file}/g" elev_parameters.txt | awk '{ if (NR > 1) print int($1); else print $0}' > ss_parameters.txt

xmin=$( awk '{ if (NR == 2) print $1}' ss_parameters.txt )
xmax=$( awk '{ if (NR == 3) print $1}' ss_parameters.txt )
ymin=$( awk '{ if (NR == 4) print $1}' ss_parameters.txt )
ymax=$( awk '{ if (NR == 5) print $1}' ss_parameters.txt )

echo "dimensions for ${reduced_ss_file}"
echo $xmin $xmax $ymin $ymax
echo


xyz2grd ${reduced_ss_file} -G${binary_ss_file}=bf -I${resolution} -R${xmin}/${xmax}/${ymin}/${ymax} -N${min_ss}

grdimage ${binary_ss_file}=bf -P -JX6i+ > quick.ps



fi


# now make the topo grid


echo
echo "creating topo file"
echo





smaller_grid_file=smaller_grid_file.grd




grdcut ${topo_location} -G${smaller_grid_file} -N-9999999 -R${long_min}/${long_max}/${lat_min}/${lat_max}



# output grid is in metres (change the -A option to -Ak for kilometers)

area_grid=area.grd

echo ${x_resolution}/${y_resolution}


grdproject ${smaller_grid_file}  -R${long_min}/${long_max}/${lat_min}/${lat_max} -Jb${center_longitude}/${center_latitude}/${southern_standard_parallel}/${northern_standard_parallel}/${map_width} -G${area_grid} -D${x_resolution}/${y_resolution} -Ae -C -V  # the -C option means that the grid is relative to the defined projection center (i.e. that is where the 0,0 coordinate is)

makecpt -Cglobe -T-10000/10000/50   > shades.cpt
grdimage -Cshades.cpt -JX6i+ area.grd -P > modern_topo.ps

# the problem with grdproject is that the dx and dy are not exactly the correct x/y spacings, probably due to a precision error. It also
# doesn't center things on the origin. I've created a program that will fix these issues

# get the x and y values. Note that this will not work in DASH shell (which is the default shell in Ubuntu and Linux Mint)

num_x=$( ncdump -h ${area_grid} | grep $'\tx = ' | sed -e 's/ //g' | sed 's/.*=//' | sed 's/;//g' )
num_y=$( ncdump -h ${area_grid} | grep $'\ty = ' | sed -e 's/ //g' | sed 's/.*=//' | sed 's/;//g' )

echo $num_x $num_y

if [ -z "$num_x" ]
then
echo "**********************************"
echo "ncdump is not extracting correctly"
echo "**********************************"

echo "area grid here?"
echo ${area_grid}
ls -l 



ncdump -h ${area_grid}
exit 0

fi
ncdump  -v x  ${area_grid}  | awk '/data:/{y=1;next}y' | sed -e 's/[x={}; ]//g' | tr -d '\n' | sed -e 's/,/\n/g' > x_values.txt

ncdump  -v y  ${area_grid}  | awk '/data:/{y=1;next}y' | sed -e 's/[y={}; ]//g' | tr -d '\n' | sed -e 's/,/\n/g' > y_values.txt




# output to an ASCII file (warning, might be huge)

xyz_file="index_topo.xyz"



# the -Cf option outputs it into array indices for reading into a Fortran program

grd2xyz ${area_grid} -V -sa -Cf > ${xyz_file}


# Since things are not in a nice grid starting at 0,0, this program converts it.

even_grid_file="even_grid.dat"

nearest_int ${num_x} ${num_y} ${resolution} ${xyz_file} ${even_grid_file}

reduced_dem_file="modern_topo.xyz"


reduce_dem ${even_grid_file} ${outline_file} ${reduced_dem_file}

binary_dem_file="modern_topo.bin"

sed -e "s/${reduced_dem_file}/${binary_dem_file}/g" elev_parameters.txt | awk  '{ if (NR > 1) print int($1); else print $0}' > elev_parameters_temp.txt

mv -f elev_parameters_temp.txt elev_parameters.txt

xmin=$( awk '{ if (NR == 2) print $1}' elev_parameters.txt )
xmax=$( awk '{ if (NR == 3) print $1}' elev_parameters.txt )
ymin=$( awk '{ if (NR == 4) print $1}' elev_parameters.txt )
ymax=$( awk '{ if (NR == 5) print $1}' elev_parameters.txt )


echo "dimensions for ${reduced_dem_file}"
echo $xmin $xmax $ymin $ymax
echo


area_grid2="area2.grd"
xyz2grd ${reduced_dem_file} -G${area_grid2} -I${resolution} -R${xmin}/${xmax}/${ymin}/${ymax} 

plot_base="base_topo2.ps"



grdimage -Cshades.cpt -JX6i+ ${area_grid2} -P >  ${plot_base}

# include GIA



gia_file="deformation.txt"

if [  -f ${gia_file} ]
then

echo "apparently the GIA file exists: " ${gia_file} 

#ls -l -t

awk -v long_min=${long_min} -v long_max=${long_max} -v lat_min=${lat_min} -v lat_max=${lat_max}  '{ if ( $1 >= long_min && $2 >= lat_min && $1 <= long_max && $2 <= lat_max) {print $1, $2, $3 }}' ${gia_file} > ${ice_model}_gia.txt 

mapproject ${ice_model}_gia.txt   -R${long_min}/${long_max}/${lat_min}/${lat_max} -F -C -Jb${center_longitude}/${center_latitude}/${southern_standard_parallel}/${northern_standard_parallel}/${map_width}  > transformed_${ice_model}_gia.txt 

blockmean transformed_${ice_model}_gia.txt  -I${resolution} -R${xmin}/${xmax}/${ymin}/${ymax} > mean_data_gia.txt

triangulate   mean_data_gia.txt -Ggia.grd -I${resolution} -R${xmin}/${xmax}/${ymin}/${ymax} 

area_grid3=gia_deformed_topo.grd

grdmath ${area_grid2} gia.grd ADD = ${area_grid3}

area_grid4=${area_grid3}

else

area_grid4=${area_grid2}


fi


binary_dem_file2="base_topo.bin"

grdreformat ${area_grid4} ${binary_dem_file2}=bf 

sed -e "s/${binary_dem_file}/${binary_dem_file2}/g" elev_parameters.txt | awk  '{ if (NR > 1) print int($1); else print $0}' > elev_parameters_temp.txt

mv -f elev_parameters_temp.txt elev_parameters.txt

# test plot

plot_topo="base_topo3.ps"



smoothing_factor=1

echo "${resolution} * ${smoothing_factor}" | bc 
 

map_width=0.5c # likely arbitrary
xtext="x"

xtickint=200000
xsubtickint=100000

ytext="y"

ytickint=200000
ysubtickint=100000


grdimage ${area_grid2} -Cshades.cpt -X5 -Y7 -JX15/15 -K -R${xmin}/${xmax}/${ymin}/${ymax}  -Ba"${xtickint}"f"${xsubtickint}":"${xtext}":/a"${ytickint}"f"${ysubtickint}":"${ytext}":WeSn -V -P  > $plot_topo


psxy ${outline_file} -K -O -R -JX -Wthickest,yellow >> ${plot_topo}
psscale -X-3 -Y-1 -D9c/-2c/15c/0.5ch -O -Ba1000:"Topography (m)": -Cshades.cpt  >> $plot_topo


