/*******************************************************************************
    * CRU CL 1.0 Data Extraction Routines  - retrieves site based readings from 0.5x0.5deg 
    * global climate data set available at http://www.cru.uea.ac.uk/cru/data/hrg.htm
    * Date   : 22/5/05
    * Copyright (C) 2005 Daniel Falster
    *
    * This program is free software; you can redistribute it and/or
    * modify it under the terms of the GNU General Public License
    * as published by the Free Software Foundation; either version 2
    * of the License, or (at your option) any later version.
    * 
    * This program is distributed in the hope that it will be useful,
    * but WITHOUT ANY WARRANTY; without even the implied warranty of
    * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    * GNU General Public License for more details.
    * 
    * The GNU General Public License can be viewed at http://www.gnu.org/copyleft/gpl.html
    * or obtained from the Free Software Foundation, 
    * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,USA.
    *
    * Contact details: Daniel Falster, dfalster@bio.mq.edu.au
    * ARC-NZ Research Network for vegetation Function
    * Dept Biological Sciences, Macquarie University 2109, Australia
    * http://www.bio.mq.edu.au/ecology/vegfunction/
    ***********************************************************************************/


/***********************************************************
INPUT FORMAT:
Climatic data files available for download from: http://ipcc-ddc.cru.uea.ac.uk/obs/get_30yr_means.htm
The files containing the observed climate data are in ASCII format with one file per year per variable. The first and second lines of the file contain information on the grid size, for example:
grd_sz xmin ymin xmax ymax n_cols n_rows n_months missing
0.5 0.25 -89.75 359.75 89.75 720 360 12 -9999
This is followed by 12 monthly grids that are n_cols by n_rows in size. Each record contains n_cols (=720) columns longitude values, format 720i5, starting at xmin (= 0.25 deg East) and ending at xmax (=0.25 deg West) The first record starts at January, ymax (=89.75 deg North) and the last record is the row for December, ymin (=-89.75 deg South). Co-ordinates represent the centre of each grid cell. Missing values (Antarctica, oceans and major inland water bodies) are assigned integer values of -9999.

Site list is in format
name    lat(-90 to 90) long(-150 to 150)     ***************/
 
#include <stdio.h>
#include <iostream>
#include <string>
#include <fstream>

using namespace std;

string read(ifstream& File, char stop);

int main(int argc, char *argv[])
{
// char ch;
 long int i,j,m, cols, rows, no_months, col_width=5;
 string line; line.reserve(4000);
 
 //open climate data file
 cout<<"Enter climate data filename: "; cin >>line;
 ifstream inFile(line.c_str());  if(!inFile) {cerr<<"File not opened"<<endl;system("pause"); exit(EXIT_FAILURE);}
 
 //open output file
 cout<<"Enter output file name: ";  cin >>line;
 ofstream outFile(line.c_str());  if(!outFile)  {cerr<<"OUTFile not opened"<<endl;system("pause"); exit(EXIT_FAILURE);}

 //open lat lon file
 cout<< "Enter lat-lon list filename: "; cin >> line;
 ifstream latlonFile(line.c_str());  if(!latlonFile) {cerr<<"File not opened"<<endl;system("pause"); exit(EXIT_FAILURE);}
 getline(latlonFile,line); //get rid of title line

 //enter climate data file description (first two lines)
 struct data_desc{string name; float val;};  data_desc file_descr[8];//data structure used for holding file description
 for(i=0; i<8; i++)  inFile >> file_descr[i].name;   
 for(i=0; i<8; i++)  inFile >> file_descr[i].val;
 read(inFile, '\n'); //read to end of line

 //print datafile description to screen
 for(i=0; i<8; i++) cout<< file_descr[i].name <<": " << file_descr[i].val<<endl;   cout<<endl;
 
//Read in climate data
 cols=(int)file_descr[5].val;  rows=(int)file_descr[6].val;    no_months=(int)file_descr[7].val;
 int* month_data =new int [no_months*rows*cols];   //3D matrix for storing climate data - 12 months/sites * (rows*coles) sites

 cout<<"\n\nPlease wait.........memory hungry task in operation"<<endl;
 for(m=0; m<no_months; m++)
   for(i=0; i< rows; i++)
     {getline(inFile,line);
      for(j=0;j<cols;j++)
            month_data [m*(cols*rows)+i*(cols)+j]=atoi((line.substr(j*col_width,col_width)).c_str());
     }

 //EXTRACT CLIMATE DATA FOR EACH SITE
 struct site_data_type{string locality; float lat; float lon; int row; int col;}; //data to hold site data
 site_data_type site;
 //Print headers to file
 outFile <<"Locality\tLat\tLon\tlat_row\tLon_col\t"; for(i=0;i<no_months; i++) outFile <<"month_"<<i+1<<"\t";   outFile <<endl;
 //Read in site details and print data to file
 while ( (latlonFile >> ws) && (! latlonFile.eof()) )
      {
      latlonFile >> site.locality >> site.lat >> site.lon;
      if(site.lat>= 90.0 || site.lat<=-90.0) {cout<<"invalid latitude value at "<<site.locality<<endl; system("pause");}
      if(site.lon>= 180.0 || site.lon<=-180.0) {cout<<"invalid longitude value at "<<site.locality<<endl; system("pause");}
      
      //calculate row & col of data from lat & long
      //formula = times by no cols per degree (i.e. 1/ grid size) * [ lattidue + half grid size - minimum centre point]
      site.row=(int)(1/file_descr[0].val*(-site.lat+file_descr[0].val/2-file_descr[2].val)); // NOTE: data grid treats north as negative, opposite to convential, so take negative of latitude
      if(site.lon>= 0.0)
            site.col=(int)(1/file_descr[0].val*(site.lon+file_descr[0].val/2-file_descr[1].val));
      else   //transform negative longitudes to range 180-360 
            site.col=(int)(1/file_descr[0].val*(360.0+site.lon+file_descr[0].val/2-file_descr[1].val));
      
      outFile << site.locality<<"\t"<<site.lat<<"\t"<<site.lon<<"\t"<<site.row<<"\t"<<site.col<<"\t";
      //find data for each site for each month & print to file
      for(m=0;m<no_months; m++)
              outFile << month_data [m*(cols*rows)+site.row*(cols)+site.col] <<"\t"; 
      outFile<<endl;
      }

 inFile.close();
 outFile.close();
 latlonFile.close();
 cout<<"Finished! Goobye.......\n"<<endl;
 
 
 return 1;
 }

//scans a sequence until stop character is reached
string read(ifstream& File, char stop)
    {string temp;
    char ch=File.get();
    while(!File.eof() && ch!=stop)  {temp+=ch; ch=File.get();}
    return temp+'\0';}