/************************************************************************************************************
* peakele -f /path/demfile [-d 0] [-r 100000]                                                               *
*                                                                                                           *
* Compiling with "cc -Wall -o  peakele  peakele.c -lgdal -lm -O2"                                           *
*                                                                                                           *
************************************************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <time.h>
#include <math.h>
#include "gdal/gdal.h"
#include "gdal/cpl_conv.h"
#include "gdal/gdal_frmts.h"

struct list_peak{
 long long int    id;
 long long int    heigherpeak_id;
 double           heigherpoint_lon,heigherpoint_lat;
 double           lon,lat,ele,isolation;
 char             elestring[100];
};    
  
   

#define PI             3.1415927
/* min isolation of a peak */
#define MINISO         100
/* difference between the heigts of DEM and peak to define a isolation */
#define MINDIFF         20


int    debuglevel=0;
time_t starttime;


void printhelp(){
 printf("Usage:   peakele -f demfile [options]\n");
 printf("         -f demfile\n         -d debuglevel (0-4, default=0)\n");
 printf("         -n maximal number of peaks (default 2000000)\n");
 printf("demfile is a geotiff (1 band int16, SRS=EPSG 4326)\n");
 printf("stdin   is a csv with \"id;lon;lat;ele\" where \n");
 printf("         ele is integer, float or empty (not mixed like \"1234m\")\n");
 printf("stdout  is a csv with \"id;lon;lat;ele;srtmele\"\n\n");
}

double get_height(double lon,double lat,double *adfGeoTransform,GDALRasterBandH hBand,long int xsize,long int ysize){
/* get height at lon/lat with elevation data from hBand (parameters in adfGeoTransform) */

 long int   xpx,ypx;
 double     h;
 int16_t    *area2x2;

 /* calculate upper left pixel next to lon/lat  */

 xpx=(long int)floor((lon-adfGeoTransform[0])/adfGeoTransform[1]);
 ypx=(long int)floor((lat-adfGeoTransform[3])/adfGeoTransform[5]);
    
/* read a 2x2 area from demfile at this pixel, so getting all 4 pixel next to lon/lat  */
/* take the max. height from these 4 values                                            */
 
 area2x2=malloc(4*sizeof(int16_t));
 GDALRasterIO(hBand,GF_Read,xpx,ypx,2,2,area2x2,2,2,GDT_Int16,0,0);
 if(debuglevel>3){printf("  Getting x/y:%ld %ld at %lf %lf h=%d %d %d %d\n",xpx,ypx,adfGeoTransform[0]+xpx*adfGeoTransform[1],adfGeoTransform[3]+ypx*adfGeoTransform[5],area2x2[0],area2x2[1],area2x2[2],area2x2[3]);}
 h=area2x2[0];
 if(area2x2[1]>h){h=area2x2[1];}
 if(area2x2[2]>h){h=area2x2[2];}
 if(area2x2[3]>h){h=area2x2[3];}
 free(area2x2);
 if(h<-32000){h=0.0;}
 return h;
}



int main(int argc, char *argv[]){

 GDALDatasetH     hDataset;
 GDALRasterBandH  hBand;
 double           adfGeoTransform[6];
 char             *demfile=NULL;
 char             *outputformat=NULL;
 char             *textrest=NULL;
 double           radius=100000.0;
 char             line[100],elestring[100];
 long long int    id;
 long int         n=0,i=0,st_s=0,st_g=0,maxnumpeaks=2000000,numpeaks=0;
 double           lon,lat,ele;
 struct list_peak *peak;

/* Parse command line parameters */

 starttime=time(NULL);
 if(argc<3){printhelp();exit(1);}
 for(i=1;i<=argc-2;i+=2){
  if     (strcmp(argv[i],"-d")==0){debuglevel=atoi(argv[i+1]);}
  else if(strcmp(argv[i],"-f")==0){demfile=argv[i+1];}          
  else if(strcmp(argv[i],"-r")==0){radius=strtod(argv[i+1],&textrest);}
  else if(strcmp(argv[i],"-n")==0){maxnumpeaks=strtod(argv[i+1],&textrest);}
  else if(strcmp(argv[i],"-o")==0){outputformat=argv[i+1];}          
  else if(strcmp(argv[i],"-h")==0){printhelp();exit(0);}          
  else if(strcmp(argv[i],"-?")==0){printhelp();exit(0);}
 }   
 if(!demfile){printhelp();exit(1);}
 if(radius<100){radius=100;}
 if(radius>500000){radius=500000;}

 peak=malloc(sizeof(struct list_peak)*maxnumpeaks); 
 if(!peak){
  fprintf(stderr,"Not enough memory for %ld peaks\n",maxnumpeaks);
  exit(1);
 }

/* open demfile, get srs, origin, pixel size and pointer to the band 1 */

 GDALRegister_GTiff();
 hDataset = GDALOpen(demfile,GA_ReadOnly);
 if( hDataset == NULL ){
  fprintf(stderr,"Cannot open File %s\n",demfile);
  exit(1);
 }
 if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None ){
  if(debuglevel>0){
   printf( "%ld Reading %s %dpx x %dpx x %d band\n",time(NULL)-starttime,demfile,GDALGetRasterXSize(hDataset),GDALGetRasterYSize(hDataset),GDALGetRasterCount(hDataset));
   printf( "Origin =(%.6f,%.6f) ",adfGeoTransform[0],adfGeoTransform[3]);
   printf( "Pixel Size=(%.6f,%.6f)\n",adfGeoTransform[1],adfGeoTransform[5]);
  }
  hBand=GDALGetRasterBand(hDataset,1);
 }else{
  fprintf(stderr,"Could not fetch transformation parameters from %s\n",demfile);
  exit(1);
 }

/* Read stdin, split each line in id,lon,lat,ele */

 while(fgets(line,100,stdin)){
  i=sscanf(line,"%lld;%lf;%lf;%[^\n]",&id,&lon,&lat,elestring);
  if(i==3){elestring[0]='\0';}
  strtok(elestring,";\n");
  ele=get_height(lon,lat,adfGeoTransform,hBand,GDALGetRasterXSize(hDataset),GDALGetRasterYSize(hDataset));
  if(i<3){id=0;}
 
/* Append peak to array of peaks */  
  
  if(id){
   if(n>maxnumpeaks){
    fprintf(stderr,"Too many peaks: %ld\n",n);
    exit(1);
   }else{
    peak[n].id=id;peak[n].lon=lon;peak[n].lat=lat;peak[n].ele=ele;peak[n].isolation=radius;peak[n].heigherpoint_lon=peak[n].heigherpoint_lat=0.0;
    strncpy(peak[n].elestring,elestring,99);
    n++;numpeaks++;
   }
  }
  st_g++;
 }
 for(n=0;n<numpeaks;n++){
  if(strstr(peak[n].elestring,";")){printf("%lld;%.7lf;%.7lf;\"%s\";%.0f\n",peak[n].id,peak[n].lon,peak[n].lat,peak[n].elestring,peak[n].ele);}
  else                             {printf("%lld;%.7lf;%.7lf;%s;%.0f\n",peak[n].id,peak[n].lon,peak[n].lat,peak[n].elestring,peak[n].ele);}
 }
 exit(0);
}