epu/project.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <unistd.h>
   4 #include <string.h>
   5 #include <epp.h>
   6 #include <eppl_ut.h>
   7 #include <projects.h>
   8 #include <getopt.h>
   9 #include <sys/mman.h>
  10 PJ *in_proj=NULL,*out_proj=NULL;
  11 int positive_west=0;/* Whether to adjust west hemisphere coords to range
  12                        180-360 from range -180-0 */
  13 EPP *infile,*outfile;
  14 /* This function does perform actual projection conversion
  15    It computes long/lat data from current cell and then computes
  16    cell on input map to get value from */
  17 int project_cell(int col, int row, int value) {
  18    UV out,ll,in;
  19    int in_row,in_col;
  20    out.u=alt_xc(outfile,col);
  21    out.v=alt_yc(outfile,row);
  22    if (out_proj) {
  23      ll=pj_inv(out,out_proj);
  24    } else {
  25        if (positive_west&&out.u>180) {
  26            out.u-=360.0;
  27        }
  28        ll.u=out.u*DEG_TO_RAD;
  29        ll.v=out.v*DEG_TO_RAD;
  30
  31    }
  32    if (ll.u == HUGE_VAL || ll.v == HUGE_VAL) {
  33      return value;
  34    }
  35    if (in_proj) {
  36    in=pj_fwd(ll,in_proj);
  37    } else {
  38        in.u=ll.u*RAD_TO_DEG;
  39        if (positive_west && in.u<0) {
  40            in.u+=360.0;
  41        }
  42        in.v=ll.v*RAD_TO_DEG;
  43    }
  44    in_row=epp_row(infile,in.v);
  45    in_col=epp_col(infile,in.u);
  46    return epp_get(infile,in_col,in_row);
  47 }
  48
  49 void help(exitcode) {
  50    printf("Usage: project [-%%] {-b basefile|-c cellsize -l xl,yb,xr,yt} [-i source_proj]\n\t"
  51    "[-p target_proj] -o outfile -O offsite {-r rows_to_cache|-x [tmpdir]} file\n"
  52   "Projection definition can be either file name or list of comma separated\n"
  53   "name=value pairs. If either input or output projection omitted, it is assumed\n"
  54   "to be geographic coords\n");
  55   exit(exitcode);
  56 }
  57 PJ* get_projection(const char *source);
  58 void expand_epp(EPP *epp,const char *tempdir,int verbose);
  59 void check_hemisphere(EPP *epp);
  60 int main(int argc,char **argv) {
  61 struct option long_options[]={
  62   {"help",0,0,1},
  63   {"version",0,0,2},
  64   {"verbose",0,0,'%'},
  65   {"base-file",1,0,'b'},
  66   {"cell-size",1,0,'c'},
  67   {"source-projection",1,0,'i'},
  68   {"output-file",1,0,'o'},
  69   {"target-projection",1,0,'p'},
  70   {"offsite",1,0,'O'},
  71   {"cache-rows",1,0,'r'},
  72   {"expand",2,0,'x'},
  73   {"limits",1,0,'l'}
  74   };
  75 char outname[1024]="project.out.epp";
  76 char tmpdir[1024]="/tmp";
  77 char *basefilename=NULL,*errptr;
  78 int limits_flag=0,result,offsite_flag=0,offsite=0,cache_size=0,expand_flag=0;
  79 double xright=0,xleft=0,ybottom=0,ytop=0,cellsize = 0;
  80 int c,index,verbose=0;
  81   while ((c=getopt_long(argc,argv,"c:i:o:p:O:r:x::l:b:%",long_options,&index))
  82              !=-1) {
  83       switch (c) {
  84           case 2: show_version("project","$Revision: 1.1 $");
  85           case '%': verbose = 1; break;
  86           case 'b' : if (cellsize !=0 || limits_flag) {
  87                          fprintf(stderr,"basefile cannot be combined with"
  88                                  "explicit coordinate parameters\n");
  89                          exit(1);
  90                       }
  91                      if (basefilename) {
  92                          fprintf(stderr,"Duplicate basefile specification\n");
  93                          exit(1);
  94                      }
  95                      basefilename=strdup(default_ext(optarg,".epp"));
  96                      break;
  97           case 'c' : if (basefilename) {
  98                        fprintf(stderr,"cannot use cell size and"
  99                                       "basefile at the same time\n");
 100                        exit(1);
 101                      }
 102                      cellsize=strtod(optarg,&errptr);
 103                      if (*errptr || cellsize<=0) {
 104                          fprintf(stderr,"invalid cellsize %s\n",optarg);
 105                          exit(1);
 106                      }
 107                      break;
 108           case  'i' : in_proj=get_projection(optarg);
 109                       if (!in_proj) exit(1);
 110                       break;
 111           case 'l' :  if (basefilename) {
 112                          fprintf(stderr,"Cannot use both basefile"
 113                                  "and limits\n");
 114                          exit(1);
 115                       }
 116                       if (limits_flag) {
 117                           fprintf(stderr,"Duplicate limit specification\n");
 118                           exit(1);
 119                       }
 120                       xleft = strtod(optarg,&errptr);
 121                       while (strchr(",;\n \t",*errptr)) errptr++;
 122                       if (!*errptr) {
 123                            fprintf(stderr,"Invalid limits");
 124                            exit(1);
 125                       }
 126                       ybottom = strtod(errptr,&errptr);
 127                       while (strchr(",;\n \t",*errptr)) errptr++;
 128                       if (!*errptr) {
 129                            fprintf(stderr,"Invalid limits");
 130                            exit(1);
 131                       }
 132                       xright = strtod(errptr,&errptr);
 133                       while (strchr(",;\n \t",*errptr)) errptr++;
 134                       if (!*errptr) {
 135                            fprintf(stderr,"Invalid limits");
 136                            exit(1);
 137                       }
 138                       ytop = strtod(errptr,&errptr);
 139                       if (*errptr) {
 140                            fprintf(stderr,"Invalid limits");
 141                            exit(1);
 142                       }
 143                       limits_flag=1;
 144                       break;
 145           case 'p'  : out_proj=get_projection(optarg);
 146                       if (!out_proj) exit(1);
 147                       if (!out_proj->inv) {
 148                           fprintf (stderr,"Cannot project - inverse "
 149                                   "transformation for output projection"
 150                                   "unavailable\n");
 151                           exit(1);
 152                       }
 153                       break;
 154          case 'o'  : strncpy(outname,default_ext(optarg,".epp"),1023);
 155                      outname[1023]=0;
 156                       break;
 157          case 'O' : offsite = strtol(optarg,&errptr,0);
 158                     if (*errptr || offsite<0 || offsite>65535) {
 159                          fprintf(stderr,"Invalid offsite value %s\n",optarg);
 160                          exit(1);
 161                     }
 162                     offsite_flag=1;
 163                     break;
 164          case 'r' : cache_size = strtol(optarg,&errptr,0);
 165                     if (*errptr || cache_size<0 || cache_size>30000) {
 166                         fprintf(stderr,"Invalid cache size:%s\n",optarg);
 167                         exit(1);
 168                     }
 169                     break;
 170          case 'x' : expand_flag=1;
 171                     if (optarg) {
 172                        strncpy(tmpdir,optarg,1023);
 173                        tmpdir[1023]=0;
 174                     }
 175                     break;
 176          default: {
 177                    help(c==1?0:1);
 178          }
 179       }
 180    }
 181    if (cache_size && expand_flag) {
 182       fprintf(stderr,"cannot use cache and expanding to temp dir"
 183               "simulateously\n");
 184       exit(1);
 185    }
 186    if (!in_proj) {
 187        if (!out_proj) {
 188            fprintf(stderr,"Neither input nor output projection specified\n");
 189            exit(1);
 190        } else {
 191         fprintf(stderr,"input projection is not specified, assuming latlong\n");
 192        }
 193    } else if (!out_proj) {
 194        fprintf(stderr,"Output projection is not specified, assuming latlong\n");
 195    }
 196    if (optind>=argc) {
 197        fprintf(stderr,"No input file specified\n");
 198        exit(1);
 199    }
 200
 201    infile=open_epp(default_ext(argv[optind],".epp"));
 202    if (!infile) {
 203         perror("open input file");
 204         exit(2);
 205    }
 206    if (!in_proj) {
 207        check_hemisphere(infile);
 208    }
 209    if (cache_size) {
 210        set_epp_cache(infile,cache_size);
 211    }
 212    if (!offsite_flag) {
 213      offsite = infile->offsite;
 214    }
 215    Create16bit=infile->kind==16;
 216    if (basefilename) {
 217        EPP *basefile;
 218        basefile=open_epp(basefilename);
 219        if (!basefile) {
 220            perror("open base file");
 221            exit(2);
 222        }
 223        outfile=creat_epp_as(outname,basefile);
 224        if (!outfile) {
 225            perror("opening output file\n");
 226            exit(1);
 227        }
 228            outfile->offsite=offsite;
 229        close_epp(basefile);
 230    } else {
 231       int rows;int cols;
 232       if (!limits_flag || cellsize==0) {
 233           fprintf(stderr,"output file geometry is not specified\n");
 234           exit(2);
 235       }
 236       rows=ceil(fabs(ytop-ybottom)/cellsize);
 237       cols=ceil(fabs(xright-xleft)/cellsize);
 238       if (ybottom<ytop) {
 239           ytop=ybottom+rows*cellsize;
 240       } else {
 241           ytop=ybottom-rows*cellsize;
 242       }
 243       if (xleft<xright) {
 244           xright=xleft+cols*cellsize;
 245       } else {
 246           xright=xleft-cols*cellsize;
 247       }
 248       outfile=creat_epp(outname,1,1,cols,rows,xleft,ytop,xright,ybottom,
 249               0,0,offsite);
 250       if (!outfile) {
 251           perror("opening output file");
 252       }
 253    }
 254    if (!out_proj) {
 255        check_hemisphere(outfile);
 256    }
 257    if (expand_flag) {
 258        expand_epp(infile,tmpdir,verbose);
 259    }
 260    install_progress_indicator(verbose?show_percent:check_int);
 261    result=clear_progress(for_each_cell(outfile,project_cell));
 262    if (result) {
 263        return abs(result);
 264    }
 265    close_epp(outfile);
 266    if (expand_flag) {
 267        infile->row=malloc(1);
 268        munmap(infile->cache,infile->width*(infile->lr-infile->fr)*
 269                sizeof(unsigned short));
 270        close(infile->cache_size);
 271    }
 272    close_epp(infile);
 273    return(0);
 274 }
 275
 276 PJ* get_projection(const char *source) {
 277     char *work_area,*token,**args;
 278     int argc,limit=16;
 279     PJ *pj;
 280     if (!strchr(source,'=')) {
 281        FILE *f=fopen(default_ext(source,".prj"),"r");
 282        int size;
 283        if (!f) {
 284           perror(default_ext(source,".prj"));
 285           exit(2);
 286        }
 287        fseek(f,0,SEEK_END);
 288        size = ftell(f);
 289        fseek(f,0,SEEK_SET);
 290        work_area=malloc(size+1);
 291        if(fread(work_area,1,size,f)!=size) {
 292            perror(default_ext(source,"prj"));
 293            exit(2);
 294        }
 295        work_area[size]=0;
 296        fclose(f);
 297     } else {
 298        work_area = strdup(source);
 299     }
 300     args=calloc(limit,sizeof(char *));
 301     argc=0;
 302     token=strtok(work_area,",;\n \r\t");
 303     do {
 304      args[argc++]=token;
 305      if (argc>=limit) {
 306         args=realloc(args,(limit+=16)*(sizeof(char *)));
 307      }
 308     } while ((token = strtok(NULL,",;\n \r\t")));
 309     pj = pj_init(argc,args);
 310     free(work_area);
 311     free(args);
 312     if (!pj) { fprintf(stderr,"Invalid projection definition: %s\n",
 313             source);
 314           exit(2);
 315     }
 316     return pj;
 317  }
 318  void position_mmap(EPP *epp,int row) {
 319      epp->row=((unsigned short*) epp->cache)+(epp->width)*(row-epp->fr);
 320      epp->currentline=row;
 321 }
 322 void expand_epp(EPP *epp,const char* tmpdir, int verbose) {
 323    int fd,i,k,n,row_size;
 324    char tempname[1048]="/tmp";
 325    if (tmpdir && *tmpdir) {
 326       strcpy(tempname,tmpdir);
 327    }
 328    strcat(tempname,"/eppXXXXXX");
 329    /* We use mkstemp rather than POSIX tmpfile here, becouse we need
 330     * to have control over place where file is created. I don't expect
 331     * that everybody have few gigabytes free in /tmp for unpacking large
 332     * map
 333     */
 334    fd=mkstemp(tempname);
 335    unlink(tempname);
 336    /*file has mode 0666, so it is better to not show it to anyone*/
 337    if (verbose) {
 338      fprintf(stderr,"unpacking...\n");
 339      install_progress_indicator(show_percent);
 340    } else {
 341      install_progress_indicator(check_int);
 342    }
 343    n=epp->lr - epp->fr;
 344    row_size=epp->width*sizeof(unsigned short);
 345    for (i=epp->fr,k=1;i<epp->lr;i++,k++) {
 346        epp_get(epp,epp->fc,i);
 347        if (write(fd,epp->row,row_size)!=row_size) {
 348            fprintf(stderr,"Not enough space for temporary file\n");
 349            exit(2);
 350        }
 351        if(EndLineProc(i,k,n)) {
 352            clear_progress(1);
 353            exit(2);
 354        }
 355    }
 356    clear_progress(0);
 357    free(epp->row);
 358    lseek(fd,0,SEEK_SET);
 359    epp->cache=mmap(NULL,row_size*n,PROT_READ,MAP_SHARED,fd,0);
 360    if (epp->cache==MAP_FAILED) {
 361        perror("mmap");
 362        exit(2);
 363    }
 364    epp->position=position_mmap;
 365    infile->cache_size=fd;
 366 }
 367
 368 void check_hemisphere(EPP *epp) {
 369     if (epp->XLeft>=0 && epp->XRight>=180) {
 370         fprintf(stderr,"Western hemisphere coords would be adjusted\n");
 371         positive_west=1;
 372     }
 373 }