mapcalc.c (8fdcbd0) - plMapcalc · Pawel Netzel

mapcalc.c

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <float.h>
#include <omp.h>
#include <gdal.h>
#include <ogr_srs_api.h>

#include "argtable3.h"
#include "libtcc.h"

#define MAX_INPUTS 256
#define MAX_OUTPUTS 64

char libtcc_path[10000] = "";

typedef int (*exec_type)(int,int,int,int,int,int,int,double*,double*, double*,double*);
typedef int (*begin_end_exec_type)(int,int,int,double*,double*);

/*
int fff(int COLS, int ROWS, int COL, int ROW, int INPNUM, int OUTNUM, int MEMNUM, double *IN, double *OUT, double *MEM, double* GEOTRANS) {
  int i;
  double v=0.0;
  for(i=0; i<10; i++) v+=sin(i);
  return 0;
}
*/
typedef struct {
  GDALDatasetH dataset_h;
  GDALRasterBandH band_h;
  int is_no_data;
  double no_data;
  int shift_cols, shift_rows;
  double *buffer;
  double *iobuffer;
} layer;

const char *opts_compress_none[] = { "TILED=NO", "BIGTIFF=YES", NULL};
const char *opts_compress_deflate[] = { "TILED=NO", "BIGTIFF=YES", "COMPRESS=DEFLATE", NULL};
const char *opts_compress_deflate2[] = { "TILED=NO", "BIGTIFF=YES", "COMPRESS=DEFLATE", "PREDICTOR=2", NULL};
const char *opts_compress_deflate3[] = { "TILED=NO", "BIGTIFF=YES", "COMPRESS=DEFLATE", "PREDICTOR=3", NULL};
const char *opts_compress_lzw[] = { "TILED=NO", "BIGTIFF=YES", "COMPRESS=LZW", NULL};
const char *opts_compress_lzw2[] = { "TILED=NO", "BIGTIFF=YES", "COMPRESS=LZW", "PREDICTOR=2", NULL };
const char *opts_compress_lzw3[] = { "TILED=NO", "BIGTIFF=YES", "COMPRESS=LZW", "PREDICTOR=3", NULL };
const char *default_type = "Float64";

layer *input_layers;
layer *output_layers;

int fl_quiet = 0;
int ninputs = 0;
int noutputs = 0;
int cols = 0;
int rows = 0;

int file_exists(const char *fname) {
  FILE *f = fopen(fname,"r");
  if(f) {
    fclose(f);
    return 1;
  }
  return 0;
}

void show_progress(int col, int cols) {
  if(fl_quiet) return;
  static int p = -1;
  int n = (int)(100.0*col/cols);
  if(p!=n) {
    p=n;
    fprintf(stderr,"\b\b\b\b%3d%%",p); fflush(stderr);
    if(p==100) fprintf(stderr,"\n");
  }
}

void show_message(FILE *f,char *message) {
  if(fl_quiet) return;
  fprintf(f,"%s\n",message);
  fflush(f);
}

void show_help(char *prog_name, void **argtable) {
  if(fl_quiet) return;

fprintf(stdout,"\nusage:\n\n%s",prog_name);
  arg_print_syntax(stdout, argtable, "\n");
  fprintf(stdout,"\n");
  arg_print_glossary_gnu(stdout, argtable);
  fprintf(stdout,"\n");
  fprintf(stdout,"   A user has to define -e or -p option!\n\n");
  fprintf(stdout,"\n");
  fprintf(stdout,"   plMapcalc version: %s\n\n", VERSION);
  fflush(stdout);
}

int set_num_threads(struct arg_int *th) {
  int i = 1;
  if(th->count>0)
    i = th->ival[0];
  omp_set_num_threads(i);
  return i;
}

int is_bbox_ok(layer *inputs, int ninputs) {
  double eps = 0.00000000001;
  double geo_transform[2][6];
  int i,j;
  
  if(ninputs<2) return 1;
  GDALGetGeoTransform(inputs[0].dataset_h,geo_transform[0]);
  for(i=1; i<ninputs; i++) {
	GDALGetGeoTransform(inputs[i].dataset_h,geo_transform[1]);
	for(j=0; j<6; j++)
	  if(fabs(geo_transform[0][j]-geo_transform[1][j])>eps)
	    return 0;
  }
  return 1;
}

int is_pojection_ok(layer *inputs, int ninputs) {
  if(ninputs<2) return TRUE;

int i = 1;
  int proj_ok = TRUE;
  char *p = (char *)GDALGetProjectionRef(inputs[0].dataset_h);
//  printf("\nProjection: %s\n\n",p);
  char *p1;
  OGRSpatialReferenceH ref = OSRNewSpatialReference(p);
  int proj = OSRIsProjected(ref);
  OGRSpatialReferenceH ref1;
  
  while(proj_ok && i<ninputs) {
	p1 = (char *)GDALGetProjectionRef(inputs[i].dataset_h);
	ref1 = OSRNewSpatialReference(p1);   
    if(proj) {
      if(!OSRIsSame(ref,ref1))
        proj_ok = FALSE;
    } else {
      if(!OSRIsSameGeogCS(ref,ref1))
        proj_ok = FALSE;
    }
    OSRDestroySpatialReference(ref1);
    i++;
  }
  
  OSRDestroySpatialReference(ref);
  return proj_ok;
}

char ** str_tokens(char *s, char sep) {
  int i, n=0;
  char *p;

char **tab = (char **)calloc(11, sizeof(char *));

for(i=0; i<11; i++) tab[i] = NULL;

n = 0;
  p = s;
  p++; // skip two starting characters. In windows, skipping drive letter and colon.
  p++;
  do {
    if(*p==sep) {
      *p = '\0';
      tab[n++] = p;
      if(n>10) {
        show_message(stderr, "Too many options per file!");
        exit(1);
      }
    }
    p++;
  }  while(*p!='\0');
  tab[n++] = p;

tab[0] = s;
  for(i=1; i<n; i++) {
    p=tab[i];
    p--;
    while(*p!='\0') 
      tab[i]=p--;
  }
  return tab;
}

char *get_file_name(const char *data) {
	int l = (int)strlen(data)+2;
	char *s = (char *)malloc(l);
	strcpy(s,data);
	char **toks = str_tokens(s,':');
    char *tok = toks[0];

char *result = NULL;
	if(tok!=NULL && *tok!='\0') {
	  l=(int)strlen(tok)+2;
	  result = (char *)malloc(l);
	  strcpy(result,tok);
	}
	free(s);
    free(toks);
	return result;
}

char *get_type_name(const char *data) {
	char *def = (char *)default_type;
	int l = (int)strlen(data)+2;
	char *s = (char *)malloc(l);
	char *result = NULL;

strcpy(s, data);
	char **toks = str_tokens(s,':');
    char *tok = toks[1];

if(tok!=NULL && strlen(tok)>0) {
	  l=(int)strlen(tok)+2;
	  result = (char *)malloc(l);
	  strcpy(result,tok);
	} else {
	  l=(int)strlen(def)+2;
	  result = (char *)malloc(l);
	  strcpy(result,def);
    }
	free(s);
    free(toks);
	return result;
}

double get_nodata_value(const char *data) {
	int l = (int)strlen(data)+2;
	char *s = (char *)malloc(l);
	strcpy(s,data);
	char **toks = str_tokens(s,':');
    char *tok = toks[2];

double result = 0.0;
	if(tok!=NULL)
	  result = atof(tok);
	free(s);
    free(toks);
	return result;
}

int get_isnodata_value(const char *data) {
	int l = (int)strlen(data)+2;
	char *s = (char *)malloc(l);
	strcpy(s,data);
	char **toks = str_tokens(s,':');
    char *tok = toks[2];
	int result = (tok!=NULL && *tok!='\0');
	free(s);
    free(toks);
	return result;
}

char *get_compression_name(const char *data) {
	int l = (int)strlen(data)+2;
	char *s = (char *)malloc(l);
	strcpy(s,data);
	char **toks = str_tokens(s,':');
    char *tok = toks[3];
	char *result = NULL;

if(tok!=NULL && strlen(tok)>0) {
	  l=(int)strlen(tok)+2;
	  result = (char *)malloc(l);
	  strcpy(result,tok);
	}
	free(s);
    free(toks);
	return result;
}

char **get_compression_opts(const char *data) {
	int l = (int)strlen(data)+2;
	char *s = (char *)malloc(l);
	strcpy(s,data);
	char **toks = str_tokens(s,':');
    char *tok = toks[3];
	char **result = NULL;
    if(tok!=NULL && *tok!='\0') {
      if(strcmp(tok,"DEFLATE")==0)
        result = (char **)opts_compress_deflate;
      else if(strcmp(tok,"LZW")==0)
        result = (char **)opts_compress_lzw;
      else if(strcmp(tok,"DEFLATE2")==0)
        result = (char **)opts_compress_deflate2;
      else if(strcmp(tok,"LZW2")==0)
        result = (char **)opts_compress_lzw2;
      else if(strcmp(tok,"DEFLATE3")==0)
        result = (char **)opts_compress_deflate3;
      else if(strcmp(tok,"LZW3")==0)
        result = (char **)opts_compress_lzw3;
      else
        result = (char **)opts_compress_none;
    }
	free(s);
    free(toks);
	return result;
}

int get_band_number(const char *data) {
	int l = (int)strlen(data)+2;
	char *s = (char *)malloc(l);
	strcpy(s,data);
	char **toks = str_tokens(s,':');
    char *tok = toks[1];
	int result = 1;
    if(tok!=NULL && *tok!='\0') 
      result = atoi(tok);
	free(s);
    free(toks);
    if(result==0) result = 1;
	return result;
}

int get_shift_cols(const char *data) {
    int l = (int)strlen(data)+2;
    char *s = (char *)malloc(l);
    strcpy(s,data);
    char **toks = str_tokens(s,':');
    char *tok = toks[2];
    int result = 0;
    if(tok!=NULL && *tok!='\0') 
      result = atoi(tok);
    free(s);
    free(toks);
    return result;
}

int get_shift_rows(const char *data) {
    int l = (int)strlen(data)+2;
    char *s = (char *)malloc(l);
    strcpy(s,data);
    char **toks = str_tokens(s,':');
    char *tok = toks[3];
    int result = 0;
    if(tok!=NULL && *tok!='\0') 
      result = atoi(tok);
    free(s);
    free(toks);
    return result;
}

layer *open_inputs(struct arg_str *input) {
  int i;
  int n = input->count;
  char *fname; 
  int band = 1;
  layer *layers = NULL;
  
  if(n>=MAX_INPUTS) {
	show_message(stderr,"Too many input files!!");
	exit(1);
  }

if(n!=0) {
    for(i=0; i<n; i++) {
      fname = get_file_name(input->sval[i]);

if(!file_exists(fname)) {
        show_message(stderr,"Input file does not exist!");
        exit(1);
      }
    }
    layers = (layer *)calloc(n,sizeof(layer));
    for(i=0; i<n; i++) {

fname = get_file_name(input->sval[i]);
      band = get_band_number(input->sval[i]);
      layers[i].shift_cols = get_shift_cols(input->sval[i]);
      layers[i].shift_rows = get_shift_rows(input->sval[i]);

layers[i].dataset_h = GDALOpenShared(fname,GA_ReadOnly);
      if(!(layers[i].dataset_h)) {
        show_message(stderr,"Problem with open input file!!");
        exit(1);
      }
      layers[i].band_h = GDALGetRasterBand(layers[i].dataset_h,band);
      if(!(layers[i].band_h)) {
        show_message(stderr,"Problem with open band from input file!!");
        exit(1);
      }

if(cols==0) {
        cols = GDALGetRasterBandXSize(layers[i].band_h);
        rows = GDALGetRasterBandYSize(layers[i].band_h);
      }

layers[i].no_data = GDALGetRasterNoDataValue(layers[i].band_h, &(layers[i].is_no_data));
      layers[i].buffer = (double *)malloc(cols*sizeof(double));
      layers[i].iobuffer = (double *)malloc(cols*sizeof(double));

fprintf(stderr,"Input layer %s, band: %d, shift: %d,%d", fname, band, layers[i].shift_cols, layers[i].shift_rows);
      if(layers[i].is_no_data==1)
        fprintf(stderr,", no_data: %lf",layers[i].no_data);
      fprintf(stderr,"\n");
    }
  }
  return layers;
}

void close_inputs(layer **layers) {
  int i;
  for(i=0; i<ninputs; i++) {
	GDALClose((*layers)[i].dataset_h);
	free((*layers)[i].buffer);
	free((*layers)[i].iobuffer);
  }
  free(*layers);
  *layers = NULL;
}

GDALDataType gdal_data_type(char *data) {
  if(data==NULL)
    return GDT_Float64;
  if(strcmp(data,"Byte")==0)
    return GDT_Byte;
  else if(strcmp(data,"UInt16")==0)
    return GDT_UInt16;
  else if(strcmp(data,"Int16")==0)
    return GDT_Int16;
  else if(strcmp(data,"UInt32")==0)
    return GDT_UInt32;
  else if(strcmp(data,"Int32")==0)
    return GDT_Int32;
  else if(strcmp(data,"Float32")==0)
    return GDT_Float32;
  else if(strcmp(data,"Float64")==0)
    return GDT_Float64;
  else
    return GDT_Float64;
}

double read_no_data(struct arg_dbl *output_no_data) {
  if(output_no_data->count==0)
    return DBL_MIN;
  return output_no_data->dval[0];
}

void create_outputs(struct arg_str *outputs, 
                   layer **out, layer *inp) {
  double geo_transform[6];
  char *proj, *type;
  int i;
  layer *o;
  
  int is_no_data;
  double no_data;
  GDALDataType data_type;
  char *fname; 
  char *cname = NULL;
  char **opts = NULL;

o = (layer *)malloc(outputs->count * sizeof(layer));

proj = (char *)GDALGetProjectionRef(inp[0].dataset_h);  
  GDALGetGeoTransform(inp[0].dataset_h,geo_transform);

for(i=0; i<outputs->count; i++) {

GDALDriverH driver = GDALGetDriverByName("GTiff");
    
    fname = get_file_name(outputs->sval[i]);
    type = get_type_name(outputs->sval[i]);
    data_type = gdal_data_type(type);
    no_data = get_nodata_value(outputs->sval[i]);
    is_no_data = get_isnodata_value(outputs->sval[i]);
    cname = get_compression_name(outputs->sval[i]);
    opts = get_compression_opts(outputs->sval[i]);

fprintf(stderr,"Output layer %s, type: %s",fname,get_type_name(outputs->sval[i]));
    if(is_no_data==1)
      fprintf(stderr,", no_data: %lf",no_data);
    if(cname!=NULL)
      fprintf(stderr,", compression: %s",cname);
    fprintf(stderr,"\n");

o[i].dataset_h = GDALCreate(driver,fname,
                                cols,rows,1,
                                data_type,
                                opts);
    free(type);
    free(fname);
    if(cname!=NULL) free(cname);
    
    if(!(o[i].dataset_h)) {
      show_message(stderr,"Problem with creating an output file!!");
      exit(1);
    }
     
    o[i].band_h = GDALGetRasterBand(o[i].dataset_h,1);

o[i].no_data = no_data;
    o[i].is_no_data = is_no_data;
    if(is_no_data)
      GDALSetRasterNoDataValue(o[i].band_h,no_data);

o[i].buffer = (double *)malloc(cols*sizeof(double));
    o[i].iobuffer = (double *)malloc(cols*sizeof(double));

GDALSetProjection(o[i].dataset_h,proj);
    GDALSetGeoTransform(o[i].dataset_h,geo_transform);
  
  }
  
  *out = o;
}

void close_outputs(layer **layers) {
  int i;
  if(noutputs==0) return;
  for(i=0; i<noutputs; i++) {
	GDALClose((*layers)[i].dataset_h);
	free((*layers)[i].buffer);
	free((*layers)[i].iobuffer);
  }
  free(*layers);
  *layers = NULL;
}
/*
void read_buffers(int ninputs, layer *layers, int row) {
  int i;
  CPLErr res;
  for(i=0; i<ninputs; i++) {

res=GDALRasterIO(layers[i].band_h, GF_Read, 0, row, cols, 1,
                       layers[i].buffer, cols, 1, GDT_Float64, 0, 0);
    if(res>2) {
      show_message(stderr, "Error in data reading!");
      exit(1);
    }
  }
}*/

void read_buffers(int ninputs, layer *layers, int row) {
  int i;
  CPLErr res;
  for(i=0; i<ninputs; i++) {

double d = 0.0;
    double *p;
    row -= layers[i].shift_rows;
    if(row<0 || row>=rows) {
      if(layers[i].is_no_data) d = layers[i].no_data;
      p = layers[i].buffer;
      for(i=0; i<cols; i++)
        *(p++) = d;
    } else {
      if(layers[i].shift_cols==0)
        res=GDALRasterIO(layers[i].band_h, GF_Read, 0, row, cols, 1,
                         layers[i].buffer, cols, 1, GDT_Float64, 0, 0);
      else if(layers[i].shift_cols<0) {
        res=GDALRasterIO(layers[i].band_h, GF_Read,0, row, cols+layers[i].shift_cols, 1,
                         layers[i].buffer-layers[i].shift_cols, cols + layers[i].shift_cols, 1, 
                         GDT_Float64, 0, 0);
        if(layers[i].is_no_data) d = layers[i].no_data;
        p = layers[i].buffer;
        for(i=0; i<-layers[i].shift_cols; i++)
          *(p++) = d;
      } else { // >0
        res=GDALRasterIO(layers[i].band_h, GF_Read, layers[i].shift_cols, row, cols-layers[i].shift_cols, 1,
                         layers[i].buffer, cols - layers[i].shift_cols, 1, GDT_Float64, 0, 0);
        if(layers[i].is_no_data) d = layers[i].no_data;
        p = layers[i].buffer + (cols - layers[i].shift_cols);
        for(i=0; i<layers[i].shift_cols; i++)
          *(p++) = d;
      }
    }
    if(res>2) {
      show_message(stderr, "Error in data reading!");
      exit(1);
    }
    printf("layer[%d], row%d; ",i,row);
  }
  printf("\n");
}

void write_buffers(int noutputs, layer *layers, int row) {
  int i;
  CPLErr res;
  if(noutputs==0) return;
  for(i=0; i<noutputs; i++) {

res=GDALRasterIO(layers[i].band_h, GF_Write, 0, row, cols, 1,
                     layers[i].buffer, cols, 1, GDT_Float64, 0, 0);
    if(res>2) {
      show_message(stderr, "Error in data writting!");
      exit(1);
    }
  }
}
/*
void read_iobuffer(layer *l, int cols, int row) {
  CPLErr res;
  res=GDALRasterIO(l->band_h, GF_Read, 0, row, cols, 1,
                   l->iobuffer, cols, 1, GDT_Float64, 0, 0);
  if(res>2) {
    show_message(stderr, "Error in data reading!");
    exit(1);
  }
}*/

void read_iobuffer(layer *l, int cols, int row) {

int i;
    CPLErr res = 0;
    double d = 0.0;
    double *p;
    row += l->shift_rows;
    if(row<0 || row>=rows) {
      if(l->is_no_data) d = l->no_data;
      p = l->iobuffer;
      for(i=0; i<cols; i++)
        *(p++) = d;
    } else {
      if(l->shift_cols==0)
        res=GDALRasterIO(l->band_h, GF_Read, 0, row, cols, 1,
                         l->iobuffer, cols, 1, GDT_Float64, 0, 0);
      else if(l->shift_cols<0) {
        res=GDALRasterIO(l->band_h, GF_Read,0, row, cols+l->shift_cols, 1,
                         l->iobuffer-l->shift_cols, cols + l->shift_cols, 1, 
                         GDT_Float64, 0, 0);
        if(l->is_no_data) d = l->no_data;
        p = l->iobuffer;
        for(i=0; i<-l->shift_cols; i++)
          *(p++) = d;
      } else { // >0
        res=GDALRasterIO(l->band_h, GF_Read, l->shift_cols, row, cols-l->shift_cols, 1,
                         l->iobuffer, cols - l->shift_cols, 1, GDT_Float64, 0, 0);
        if(l->is_no_data) d = l->no_data;
        p = l->iobuffer + (cols - l->shift_cols);
        for(i=0; i<l->shift_cols; i++)
          *(p++) = d;
      }
    }
    if(res>2) {
      show_message(stderr, "Error in data reading!");
      exit(1);
    }
}

void write_iobuffer(layer *l, int cols, int row) {
  CPLErr res;
  res=GDALRasterIO(l->band_h, GF_Write,0, row, cols, 1,
              	   l->iobuffer, cols, 1, GDT_Float64, 0, 0);
  if(res>2) {
    show_message(stderr, "Error in data writting!");
    exit(1);
  }
}

void swap_buffers(layer *l) {
  double *d;
  d = l->buffer;
  l->buffer = l->iobuffer;
  l->iobuffer = d;
}

char *empty = "\0";

void get_code(struct arg_file *file, struct arg_str *prog, int *is_file, long *len, char **buffer) {
  *len = 0;
  *is_file = 0;
  
  if(file->count==0 && prog->count==0)
    return;
  if(file->count>0) {
    if(strcmp("mc",file->extension[0])==0)
      return;
    FILE *f = fopen(file->filename[0],"r");
    if(f) {
      fseek(f,0,SEEK_END);
      *len=ftell(f);
      fseek(f,0,SEEK_SET);
      *buffer=(char*)malloc(*len);
      *is_file=1;
      fread(*buffer,1,*len,f);
      fclose(f);
    } else {
      show_message(stderr, "Cannot read a macro file!");
      exit(2);
    }
  } else {
    *len=(long)strlen(prog->sval[0]);
    *buffer=(char *)prog->sval[0];
  }
}

char *prepare_code(struct arg_file *file, struct arg_str *prog,
                   struct arg_file *file_begin, struct arg_str *prog_begin,
                   struct arg_file *file_end, struct arg_str *prog_end) {
  char prog_template[] = 
       "#include <math.h>\n" 
       "#include <stdio.h>\n" 
       "int start_count = 1;\n"
       "int restart_flag = 0;\n"
       "void RESTART() {start_count++; restart_flag=-1;}\n" 
       "int ITERATION() {return start_count;}\n" 
       "int exec_begin(int COLS, int ROWS, int MEMNUM, double *MEM, double* GEOTRANS) {\n"
       "  restart_flag = 0;\n"
       "  %s\n"
       "  return 0;\n"
       "}\n"
       "int exec_cell(int COLS, int ROWS, int COL, int ROW, int INPNUM, int OUTNUM, int MEMNUM, double *IN, double *OUT, double *MEM, double* GEOTRANS) {"
       "  restart_flag = 0;\n"
       "  %s\n"
       "  return restart_flag;\n"
       "}\n"
       "int exec_end(int COLS, int ROWS, int MEMNUM, double *MEM, double* GEOTRANS) {\n"
       "  restart_flag = 0;\n"
       "  %s\n"
       "  return restart_flag;\n"
       "}\n\0";
  char *code = NULL;
  char *code_cell = empty;
  char *code_begin = empty;
  char *code_end = empty;
  long l = (long)strlen(prog_template);
  long ll;

int is_file_cell = 0;
  int is_file_begin = 0;
  int is_file_end = 0;

get_code(file_begin, prog_begin, &is_file_begin, &ll, &code_begin);
  l+=ll;
  get_code(file, prog, &is_file_cell, &ll, &code_cell);
  l+=ll;
  get_code(file_end, prog_end, &is_file_end, &ll, &code_end);
  l+=ll;
  l++;
  code=malloc(l);
  sprintf(code, prog_template, code_begin, code_cell, code_end);
  if(is_file_begin)
    free(code_begin);
  if(is_file_cell)
    free(code_cell);
  if(is_file_end)
    free(code_end);

return code;
}

void compile_code(TCCState* tcc, char *code) {

if (!tcc) {
    show_message(stderr, "Cannot create executable code (stage 1)");
    exit(1);
  }

char *path = malloc(strlen(libtcc_path) + 30);

tcc_set_lib_path(tcc, libtcc_path);

strcpy(path, libtcc_path);
  strcat(path, "include");
  tcc_add_include_path(tcc, path);

strcpy(path, libtcc_path);
  strcat(path, "lib");
  tcc_add_library_path(tcc, path);

free(path);
// #else
  /* Linux 64 section */
  //tcc_add_library_path(tcc, "/usr/share/mapcalc/");
//  tcc_set_lib_path(tcc, "/usr/local/share/plmapcalc");
//  tcc_add_include_path(tcc, "/usr/loacl/share/plmapcalc/include");
//  tcc_add_include_path(tcc, "/usr/local/share/plmapcalc/tcc_include");
//  tcc_set_lib_path(tcc, "/usr/local/lib/tcc");
//  tcc_add_include_path(tcc, "/usr/local/lib/tcc/include");
//  tcc_add_include_path(tcc, "/usr/local/include");
// #endif

tcc_set_output_type(tcc, TCC_OUTPUT_MEMORY);

if (tcc_compile_string(tcc, code) == -1) {
    show_message(stderr, "Cannot create executable code (stage 2)");
    exit(1);
  }

if (tcc_relocate(tcc, TCC_RELOCATE_AUTO) < 0) {
    show_message(stderr, "Cannot create executable code (stage 3)");
    exit(1);
  }
}

exec_type get_cell_function(TCCState* tcc) {

exec_type e = tcc_get_symbol(tcc, "exec_cell");
  if (!e) {
    show_message(stderr, "Cannot create CELL executable");
    exit(1);
  }
  
  return e;
}

begin_end_exec_type get_begin_function(TCCState* tcc) {

begin_end_exec_type e = tcc_get_symbol(tcc, "exec_begin");
  if (!e) {
    show_message(stderr, "Cannot create BEGIN executable");
    exit(1);
  }
  
  return e;
}

begin_end_exec_type get_end_function(TCCState* tcc) {

begin_end_exec_type e = tcc_get_symbol(tcc, "exec_end");
  if (!e) {
    show_message(stderr, "Cannot create END executable");
    exit(1);
  }
  
  return e;
}

int is_null(layer *l, double v) {
  if((l->is_no_data) && (v==l->no_data))
    return TRUE;
  return FALSE;
}

void save_memory_cells(char *fname, double *mem, int nmem, int withzeros) {
  int i;
  FILE *f = fopen(fname,"w");
  if(f==NULL) {
    fprintf(stderr, "\022[1;31mWarning: Can not store memory data into file: %s\n",fname);
  } else {
    for(i=0; i<nmem; i++)
      if(mem[i]==0.0) {
        if(withzeros!=0)
          fprintf(f, "%d %.18lf\n", i, mem[i]);
      } else
        fprintf(f, "%d %.18lf\n", i, mem[i]);
    fclose(f);
  }
}

void read_memory_cells(char *fname, double *mem, int nmem) {
  int i,k;
  double d;
  FILE *f = fopen(fname,"r");
  if(f==NULL) {
    fprintf(stderr, "\022[1;31mWarning:\033[0m Can not read memory data from file: %s\n",fname);
  } else {
    i=0;
    while(!feof(f)) {
      k = fscanf(f, "%d %lf\n", &i, &d);
      if(k==2 && i>=0 && i<nmem)
        mem[i]=d;
    }
    fclose(f);
  }
}

int read_max_memory(char *fname) {
  int i,k, max;
  double d;
  FILE *f = fopen(fname,"r");
  i=0;
  max = -1;
  while(!feof(f)) {
    k = fscanf(f, "%d %lf\n", &i, &d);
    if(k==2 && i>max) max = i;
  }
  fclose(f);
  return max+1;
}

int main(int argc,char **argv) {
    TCCState* tcc; 
    exec_type exec;
    begin_end_exec_type exec_begin = NULL;
    begin_end_exec_type exec_end = NULL;
    char *code;
    int nerr;
    int row;
    double *in_data;
    double *out_data = NULL;
    double *memory = NULL;
    int memnum = 0;
    int ip;
    int remove_nan = FALSE;

/* args */
    struct arg_lit  *help = arg_lit0("h", "help", "program usage");
    struct arg_lit  *quiet = arg_lit0("q", "quiet", "quiet mode");
    struct arg_lit  *check = arg_lit0(NULL, "check", "display parameters and check their correctness");
    struct arg_lit  *overwrite = arg_lit0("f", "force", "force to overwrite output file");
    struct arg_lit  *nan = arg_lit0(NULL, "use-nan", "force to treat NAN as a value (default: no-data)");
    struct arg_int  *threads = arg_int0("t", "threads", "<n>","number of threads");
    struct arg_str  *input = arg_strn("i", "input", "file.tif[[[:band]:shit_x]:shif_y]",1,MAX_INPUTS,"input layer(s) (GeoTIFF or any other GDAL raster format), optionally followed by the band number, horizontal shift by the specified number of cells (a negative value means a rightward shift), vertical shift by the specified number of cells (a negative value means a downward shift)");
    struct arg_str  *output = arg_strn("o", "output", "file.tif[:type[:no_data:[compression]]]",0,MAX_OUTPUTS,"output layer(s): file name, data type (Byte, Int16, UInt16, Int32, UInt32, Float32, default Float64), no_data value (default 0.0), compress (NONE, DEFLATE, DEFLATE2, DEFLATE3, LZW, LZW2, LZW3). (GeoTIFF)");
    struct arg_int  *mem = arg_int0("m", "memory", "<n>","number of memory cells to store temporary data");
    struct arg_str  *memstore = arg_str0("s", "memory-store", "<file name>","file to store memory cells (TXT)");
    struct arg_str  *memread  = arg_str0("r", "memory-read", "<file name>","file to read memory cells (TXT)");
    struct arg_lit  *memzero  = arg_lit0("0", "store-zeroes", "store all memory cells including zeros (defalt: no)");
    struct arg_str  *prog = arg_str0("e", "execute", "' code '","code to execute");
    struct arg_file *file = arg_file0("p", "program", "<file>.mc","file with code to execute");
    struct arg_str  *prog_begin = arg_str0(NULL, "execute-begin", "' code '","code to execute at the start");
    struct arg_file *file_begin = arg_file0(NULL, "program-begin", "<file>.mc","file with code to execute at the start");
    struct arg_str  *prog_end = arg_str0(NULL, "execute-end", "' code '","code to execute at the end");
    struct arg_file *file_end = arg_file0(NULL, "program-end", "<file>.mc","file with code to execute at the end");
    struct arg_str  *path = arg_str0(NULL, "path", "<directory path>","path to TCC library and resurces");
    struct arg_end  *end = arg_end(20);
    void *argtable[] = {help, quiet, check, overwrite, nan,
		                threads,
		                input, output, 
		                mem, memread, memstore, memzero,
		                prog, file,
		                prog_begin, file_begin,
		                prog_end, file_end,
		                path,
		                end};

nerr = arg_parse(argc,argv,argtable);

if(nerr!=0 || help->count>0) {
      show_help(argv[0],argtable);
      arg_freetable(argtable,sizeof(argtable)/sizeof(argtable[0]));
      return 1;
    }

if(prog->count==0 && file->count==0 &&
       prog_begin->count==0 && file_begin->count==0 &&
       prog_end->count==0 && file_end->count==0) {
      printf("\nNo program code to execute!\n\n");
      arg_freetable(argtable,sizeof(argtable)/sizeof(argtable[0]));
      return 1;
    }

if(path->count>0) {
      strcpy(libtcc_path, path->sval[0]);
    } else {
      char *p = getenv("PLMAPCALCPATH");
      if(p!=NULL && strcmp(p,"")!=0)
        strcpy(libtcc_path, p);
    }

remove_nan = (nan->count==0);

fl_quiet = (quiet->count>0);
      
      ninputs = input->count;
      noutputs = output->count;

in_data = (double *)malloc(ninputs*sizeof(double));
      if(noutputs>0) {
        out_data = (double *)malloc(noutputs*sizeof(double));

/* Checking of the existance of output files */
        if(overwrite->count==0 && file_exists(output->sval[0])) {
          for(ip=0; ip<noutputs; ip++) {
            char *s = get_file_name(output->sval[ip]);
            if(file_exists(s)) {
              show_message(stderr,"Output file already exists. Use -f flag to force owerwrite.");
              exit(1);
            }
            free(s);
          }
        }
      }

/* Compile code to run */
      code = prepare_code(file, prog, file_begin, prog_begin, file_end, prog_end);
      if(code==NULL) {
        show_message(stderr, "Could not read code to execute.");
        exit(1);
      }

tcc = tcc_new();

compile_code(tcc,code);
      exec_begin = get_begin_function(tcc);
      exec = get_cell_function(tcc);
      exec_end = get_end_function(tcc);

/* Open input layers, checking projection and extent */
      GDALAllRegister();
      input_layers = open_inputs(input);
      if(!is_pojection_ok(input_layers,ninputs)) {
        show_message(stderr, "All input files should have the same projection.");
        exit(1);
      };
      if(!is_bbox_ok(input_layers,ninputs)) {
        show_message(stderr, "All input files should have the same extent and resolution.");
        exit(1);
      };

/* Creating output layers */
      if(noutputs>0) {
        create_outputs(output,&output_layers,&(input_layers[0]));
      }

if(mem->count>0) {
        memnum = mem->ival[0];
        if(memnum<1) 
          memnum=0;
        else
          memory = (double *)calloc(memnum,sizeof(double));
      }

if(memread->count>0) {
        if(file_exists((char *)(memread->sval[0]))) {
          fprintf(stderr, "Memory data will be read from file: %s\n",memread->sval[0]);
          if(memory==NULL) {
		    memnum = read_max_memory((char *)(memread->sval[0]));
            if(memnum==0) {
              show_message(stderr, "File with memory data is not properly formated!");
              exit(1);
            } else 
              memory = (double *)calloc(memnum,sizeof(double));
          }
          read_memory_cells((char *)(memread->sval[0]), memory, memnum);
        } else {
          show_message(stderr, "File with memory data does not exist!");
          exit(1);
        }
      }

if(memory!=NULL && memstore->count>0)
        fprintf(stderr, "Memory data will be stored into file: %s\n",memstore->sval[0]);

if(threads->count>0 && threads->ival[0]>0 && threads->ival[0]<100) {
        if(threads->ival[0]>ninputs+noutputs+1 && memory!=NULL)
		  omp_set_num_threads(ninputs+noutputs+1);
        else
		  omp_set_num_threads(threads->ival[0]);
      } else 
        omp_set_num_threads(1);

int data_reset;

double *geotransform = malloc(6*sizeof(double));
	  GDALGetGeoTransform(input_layers[0].dataset_h, geotransform);

double t1 = omp_get_wtime();

data_reset = exec_begin(cols,rows,memnum,memory,geotransform);
      do {
        data_reset = 0;
        show_progress(0,rows);

int i;
        for(i=0; i<ninputs; i++) {
		  read_iobuffer(&(input_layers[i]), cols, 0);
          swap_buffers(&(input_layers[i]));
        }

for(row=0; row<rows; row++) {
          show_progress(row,rows);

#ifndef WIN64
		  /* Linux 64 section - OpenMP 5*/
#pragma omp parallel default(none) shared(ninputs, noutputs, row, rows, cols, memory, memnum, input_layers, output_layers, data_reset, remove_nan, exec, geotransform) 
{
#pragma omp single 
{
#pragma omp task shared(ninputs, noutputs, row, rows, cols, memory, memnum, input_layers, output_layers, data_reset, remove_nan, exec, geotransform) 
{
#endif

#pragma omp parallel for default(shared) schedule(dynamic,100) if(memory==NULL)
          for(int col=0; col<cols; col++) {

if(data_reset==0) {
              double in[MAX_INPUTS];
              double out[MAX_OUTPUTS];
              double v;
              layer *l = NULL;
              int data_ok = TRUE;
              int i, res;

for(i=0; data_ok && i<ninputs; i++) { 
                l = input_layers + i;
                v = l->buffer[col];
                if(is_null(l,v) || (remove_nan && (isnan(v) || isinf(v)))) {
                  data_ok = FALSE;
                  break;
                } else
                  in[i] = v;
              }

if(data_ok) {
                if(noutputs>0)
                  memset(out,0,noutputs*sizeof(double));

res = exec(cols,rows,col,row,ninputs,noutputs,memnum,in,out,memory,geotransform);
                if(data_reset==0 && res!=0) {
                  #pragma omp critical 
                  {
                    data_reset = res;
                  }
                }
                for(i=0; i<noutputs; i++)
                  output_layers[i].buffer[col] = out[i];
              } else {
                for(i=0; i<noutputs; i++) {
                  l = output_layers + i;
                  l->buffer[col] = l->no_data;
                }
              }
            }
          }

#ifndef WIN64
}
#endif
          int i;
          if(row<rows-1)
            for(i=0; i<ninputs; i++)
#ifndef WIN64
#pragma omp task shared(input_layers,i,cols,row)
#endif
{
			  read_iobuffer(&(input_layers[i]), cols, row+1);
}
          if(row>0)
			for(i=0; i<noutputs; i++)
#ifndef WIN64
#pragma omp task shared(input_layers,i,cols,row)
#endif
{
			  write_iobuffer(&(output_layers[i]), cols, row-1);
}
#ifndef WIN64
#pragma omp taskwait
}
}
#endif
	  for (i = 0; i<ninputs; i++)
  		  swap_buffers(&(input_layers[i]));
        for(i=0; i<noutputs; i++)
  		  swap_buffers(&(output_layers[i]));
      }
      for(i=0; i<noutputs; i++)
		write_iobuffer(&(output_layers[i]), cols, rows-1);

show_progress(rows,rows);
      if(data_reset==0)
        data_reset = exec_end(cols,rows,memnum,memory,geotransform);

} while(data_reset!=0);

if(memory!=NULL) {
        if(memstore->count>0)
          save_memory_cells((char *)(memstore->sval[0]), memory, memnum, memzero->count);
        free(memory);
      }

close_inputs(&input_layers); 
      close_outputs(&output_layers);

double t2 = omp_get_wtime();
if(quiet->count==0) {
  fprintf(stderr,"\nTime: %lf\n\n",t2-t1);
  fprintf(stderr,"  plMapcalc version: %s\n\n", VERSION);
}
      tcc_delete(tcc);
//      free(in);
//      free(out);
      free(in_data);
      free(out_data);
      free(code);

arg_freetable(argtable,sizeof(argtable)/sizeof(argtable[0]));

return 0;
}

plMapcalc

layer shifting added