/**
 * \file asc2drainage2asc.c
 *
 * $Id$
 *
 * Taken from drainage.c which was originally written by Stefanie Rost,
 * and maintained by Werner von Bloh.
 *
 * Reads a drainage map in ASCII format, e.g. DDM30.asc,
 * and a LPJ grid.bin file.
 * Produces a new drainage map in ASCII format, adapted to the
 * land/sea boundary that is given by the grid description file.
 * That can then be use to edit the drainage map in such a way that
 * the rivers do not accidentally run out before reaching the coast.
 *
 * Encoding of the original DDM30.asc :
 * https://www.uni-frankfurt.de/45217896/3_drainage_direction_map
 * The map is generated with a cellsize of 30 arc minutes and is
 * distributed as a zipped ASCII-grid that can be easily imported in
 * most GIS-software that support rasters or grids. The direction of
 * outflow from cell is according to the ArcInfo/ArcView convention as
 * follows: "0" is assigned to internal sinks and to cells draining
 * into the ocean while "1", "2", "4", "8", "16", "32", "64" and "128"
 * standing for the outflow directions E, SE, S, SW, W, NW, N or NE.  
 * -9 : missing value / ocean cell
 *
 * Encoding of the new .asc map file:
 * -9 : missing value / ocean cell
 * -88 : missing value for land cell (according to LPJ grid.bin)
 * "1", "2", "4", "8", "16", "32", "64", "128"
 * "E", "SE","S", "SW", "W", "NW", "N",  "NE"
 * 
 */
#include "lpj.h"
#include <assert.h>

int main(int argc,char **argv)
{  
  FILE *ifp;
  FILE *mfp;
  FILE *ofp;
  Header header_grid;

  // Metadata from .asc file
  int ncols, nrows;
  float xllcorner, yllcorner, cellsize; // degrees
  int nodata;

  int i, j;  
  Bool swap;
  Intcoord coord;
  int version;

  int *inputmap, *outputmap;

  version=READ_VERSION;
  if (argc>1 && !strcmp(argv[1],"-longheader")) {
    version=2;
    argc--;
    argv++;
  } 

  /* Parse command line */
  /* e.g. "drainage re-ordered_grid DDM30.asc output" */
  if (argc!=4) {
    fprintf(stdout,"USAGE: %s [-longheader] re-ordered_grid DDM30.asc outfile\n", argv[0]);
    exit(99);
  }
 
  /* Open in- & output file */
  if ((ifp = fopen(argv[1],"rb")) == NULL) {
    fprintf(stderr,"Warning: File open failed on grid file '%s': %s\n",
            argv[1],strerror(errno));
    exit(1);
  }

  if ((ofp=fopen(argv[3],"w")) == NULL) {
    fprintf(stderr,"Warning: File open failed on output file '%s': %s\n",
            argv[3],strerror(errno));
    exit(1);
  }

  if ((mfp=fopen(argv[2],"r")) == NULL) {
    fprintf(stderr,"Warning: File open failed on drainage input '%s': %s\n",
            argv[2],strerror(errno));
    exit(1);
  } 
  /* read metadata lines from .asc file */
  fscanf(mfp,"%*s %d",&ncols);
  fscanf(mfp,"%*s %d",&nrows);
  fscanf(mfp,"%*s %f",&xllcorner);
  fscanf(mfp,"%*s %f",&yllcorner);
  fscanf(mfp,"%*s %f",&cellsize);
  fscanf(mfp,"%*s %d",&nodata);

  inputmap = newvec(int, ncols*nrows);
  for (j = 0 ; j < nrows; j++) {
    for (i = 0; i < ncols; i++) {
      if (fscanf(mfp, "%d", inputmap+i+j*ncols) != 1) {
        fprintf(stderr, "cannot read code for cell %d,%d from %s\n", i, j, argv[2]);
        perror("fscanf failed");
        exit(1);
      }
    }
  }
  outputmap = newvec(int, ncols*nrows);
  for (i = 0; i < ncols*nrows; i++) outputmap[i] = nodata;
  fprintf(ofp, "ncols         %d\n", ncols);
  fprintf(ofp, "nrows         %d\n", nrows);
  fprintf(ofp, "xllcorner     %f\n", xllcorner);
  fprintf(ofp, "yllcorner     %f\n", yllcorner);
  fprintf(ofp, "cellsize      %f\n", cellsize);
  fprintf(ofp, "NODATA_value  %d\n", nodata);

  if (freadheader(ifp,&header_grid,&swap,LPJGRID_HEADER,&version)) {
    fclose(ifp);
    fail(23,FALSE,"Invalid header in re-ordered grid files ");
  }

  for(i=0; i<header_grid.ncell; i++) {
    int ilon, jlat;
    readintcoord(ifp, &coord, swap);
    /* printf("%d %d\n",lpjlon[i],lpjlat[i]); */
    ilon = (int)((coord.lon / 100.0 - xllcorner) / cellsize);
    jlat = nrows -1 - (int)((coord.lat / 100.0 - yllcorner) / cellsize); // invert lat axis
    assert(0 <= ilon); assert(ilon < ncols);
    assert(0 <= jlat); assert(jlat < nrows);
    outputmap[ilon+jlat*ncols] = -88; // missing value for land cell (according to LPJ grid.bin)
    //printf("%6.2f %6.2f = %3d %3d <- -8\n", coord.lon / 100.0, coord.lat / 100.0, ilon, jlat);
  }
  fclose(ifp);

  // copy all non-missing-values from input map
  for (j = 0 ; j < nrows; j++) {
    for (i = 0; i < ncols; i++) {
      if (inputmap[i+j*ncols] >= 0) 
        outputmap[i+j*ncols] = (inputmap[i+j*ncols] >= 0);
    }
  }

  for (j = 0 ; j < nrows; j++) {
    for (i = 0; i < ncols; i++) {
      fprintf(ofp, "%4d ", (outputmap[i+j*ncols]));
    }
    fprintf(ofp, "\n");
  }

  fclose(ofp);

  return EXIT_SUCCESS;
}