# Set the PROJ_LIB environment variable
 ### NEEDED ? export PROJ_LIB=/opt/anaconda3/envs/gdal/share/proj

# 738  bash interpolate.sh --lat=16 --lon=-90 --res=0.2 --var=3 --range=4 --power=4 --smooth=0
# 739  gdalwarp /var/www/nube3.on.gt/data/test/interpolated.tif -r cubicspline -tr 0.015 0.015 -dstnodata -99 /var/www/nube3.on.gt/data/test/interpolatedr.tif -overwrite 
# 740  gdalwarp -cutline /vsicurl/https://raw.githubusercontent.com/Countly/countly-server/44e2c20534dac663e2d3be13d9bd7eb34fef10b2/frontend/express/public/geodata/region/GT.geojson /var/www/nube3.on.gt/data/test/interpolatedr.tif /var/www/nube3.on.gt/data/test/gt.tif -overwrite

echo $0
cd "$(dirname "$0")"
echo $(date)

. ./config.cfg 

mkdir -p $web_path ##create folder if not exists

# Define interpolation parameters
while [ $# -gt 0 ]; do
    case "$1" in
        -lat=*|--lat=*) lat="${1#*=}" ;;
        -lon=*|--lon=*) lon="${1#*=}" ;;
        -res=*|--res=*) res="${1#*=}" ;;
        -var=*|--var=*) var="${1#*=}" ;;
        -points_url=*|--points_url=*) points_url="${1#*=}" ;;
        -filename=*|--filename=*) filename="${1#*=}" ;;
        -cut=*|--cut=*) cut="${1#*=}" ;;
        -range=*|--range=*) range="${1#*=}" ;;
        -power=*|--power=*) power="${1#*=}" ;;
        -smooth=*|--smooth=*) smooth="${1#*=}" ;;
        *) echo "Unknown parameter: $1" ;;
    esac
    shift
done

# Use the parsed parameters
echo "Latitude --lat= $lat"
echo "Longitude --lon= $lon"
echo "Resolution --res= $res"
echo "Variable_id --var= $var"
echo "Points URL  (must be geojson) --points_url=$points_url"
echo "filename --filename= $filename"
echo "cut geojson --cut= $cut"
echo "range in deg --range= $range"
echo "idw power --power= $power"
echo "idw smoothing --smooth= $smooth"

# Calculate values using bc
xmin=$(echo "$lon - $range" | bc)
xmax=$(echo "$lon + $range" | bc)
ymin=$(echo "$lat - $range" | bc)
ymax=$(echo "$lat + $range" | bc)


algorithm="invdist:power=$power:smothing=$smooth:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0"       # Interpolation algorithm (inverse distance squared)

if [ "$var" -eq 7 ]; then
    algorithm="nearest:radius=1:nodata=-99.9"
fi

# Use the algorithm variable here
echo "The algorithm is: $algorithm"

output="$filename.geojson"             # Output shapefile
output_raster="$filename"    # Output interpolated raster

echo "Points: $output. Output Raster: $output_raster"

rm $web_path/$output

# Check if the points URL is set and has a length greater than 5
if [[ -n "$points_url" && ${#points_url} -gt 5 ]]; then
        # Download the points from the URL using curl
        curl -o "$web_path/$output" "$points_url"

        # Check if the download was successful (exit code 0)
        if [[ $? -eq 0 ]]; then
            echo "Downloaded points from $points_url to $output.  'var' parameter will be ignored."
        else
            echo "Error downloading points from $points_url"
            exit 1
        fi
else
        # Use ogr2ogr to query data from the database and create a shapefile
        # Define your database connection parameters
        db_host="154.12.253.150"
        db_port="11014"
        db_name="climahub"
        db_user="climahub"
        db_password="climahub\$2021"

        $gdal_path/ogr2ogr -f "GEOJSON" $web_path/$output PG:"host=$db_host port=$db_port dbname=$db_name user=$db_user password=$db_password" \
        -sql "SELECT station_id,station_name,ST_SetSRID(ST_MakePoint(station.longitude, station.latitude), 4326) AS geometry, \
        data_value::numeric valor FROM data_last, station \
        WHERE data_date >= now() - interval '30 MINUTES' AND variable_id = $var AND station_id = station.id AND 
        ST_DWithin(
            ST_MakePoint($lon,$lat)::geography,
            ST_SetSRID(ST_MakePoint(station.longitude, station.latitude), 4326)::geography,
            $range*100000 --metros
        )"
fi  

# Use gdal_grid for interpolation
set -x
## a ver si se quita el sql statement $gdal_path/gdal_grid -l sql_statement -zfield valor -a_srs EPSG:4326 -of GTiff -ot Float32 -a $algorithm -txe $xmin $xmax -tye $ymax $ymin -tr $res $res -ot Float32 $web_path/$output $web_path/$output_raster.tif
$gdal_path/gdal_grid -zfield valor -a_srs EPSG:4326 -of GTiff -ot Float32 -a $algorithm -txe $xmin $xmax -tye $ymax $ymin -tr $res $res -ot Float32 $web_path/$output $web_path/$output_raster.tif
$gdal_path/gdalwarp $web_path/$output_raster.tif -r cubicspline -tr 0.015 0.015 -dstnodata -99 $web_path/$output_raster\_r.tif -overwrite 
$gdal_path/gdalwarp -cutline /vsicurl/$cut $web_path/$output_raster\_r.tif $web_path/$output_raster\_rc.tif -overwrite
$gdal_path/gdalwarp -dstnodata -99 -cutline /vsicurl/$cut $web_path/$output_raster.tif $web_path/$output_raster\_c.tif -overwrite

##png colors rainfall
colorfile="colorsGnYlRd.txt"

# Temperatures
if [[ "$var" -eq 3 || "$var" -eq 4 ]] || [[ $filename =~ ^(.*)(TMIN|TMAX)(.*)$ ]]; then
    colorfile="colorsBuYlRd.txt"  # Use for TMIN or TMAX or specific var values
fi

# Humidity
if [[ "$var" -eq 15 ]]; then
    colorfile="colorsOrGn.txt"  # Use for specific var value
fi

# Solar rad
if [[ "$var" -eq 16 ]]; then
    colorfile="colorsGyYlRd.txt"  # Use for specific var value
fi


# Wind
if [[ "$var" -eq 8 || "$var" -eq 9 ]]; then
    colorfile="colorsBuGnRdPu.txt"  # Use for specific var values
fi


$gdal_path/gdalwarp -t_srs EPSG:3857 -r cubicspline $web_path/$output_raster\_rc.tif $web_path/$output_raster\_rc3857.tif -overwrite
$gdal_path/gdaldem color-relief $web_path/$output_raster\_rc3857.tif -nearest_color_entry $colorfile $web_path/$output_raster.png -alpha

label="Gen. $(date +'%Y-%m-%d_%H:%Mz') $output_raster" 
convert $web_path/$output_raster.png -gravity South -font arial -pointsize 10 -fill white -annotate +0+50 "$label" $web_path/$output_raster.png


set +x