dll/fgisRaster.c

   1 /* This is an implementation for EPP handling in fgis*/
   2 #include <stdlib.h>
   3 #include <tcl.h>
   4 #include <epp.h>
   5 #include <string.h>
   6 #include <strings.h>
   7 #include <reclass.h>
   8 #include <math.h>
   9 #include <eppl_ut.h>
  10 #include <sys/stat.h>
  11 #include <unistd.h>
  12 #include <errno.h>
  13 #include "fgis.h"
  14 #include "fgisInt.h"
  15 #include "fgisEppEdit.h"
  16 /*
  17  * Declaration
  18  *
  19  */
  20
  21
  22 int Fgis_CreateNewRaster(ClientData data,Tcl_Interp *interp,
  23      int argc,char **argv);
  24
  25
  26 typedef int (Fgis_CmdProc)(RASTER_OBJECT handle,Tcl_Interp *interp,
  27                           int argc, char **argv);
  28 Fgis_CmdProc  SaveRaster, ReturnBaseRaster,
  29       ReturnCellValue, ReturnLimits, RasterComments, ChangeReclass,
  30       ReturnRow, ReturnCol, DumpReclass, ReturnBPP, OffsiteValue,
  31       XLeft, XRight, YTop, YBottom,
  32       ShiftCoords, ListClasses, CountClass, FindUnused, CalcExtents,
  33       SetRasterCache, GetMinValue, GetMaxValue, CellUnit,
  34       ModifyCell, PlotLine, PlotFrame, PlotBox, PlotContour, FillContour,
  35       PlotCircle, ReturnCell, ReturnCellLimits,
  36       CountContour, CountTranssect;
  37
  38 RECLASS Fgis_ListToReclass(Tcl_Interp *interp,char *list,RECLASS src,int size);
  39 RECLASS Fgis_ParseReclass(Tcl_Interp *interp,char *program,RECLASS src,int size);
  40
  41
  42 int DumpReclassToList(Tcl_Interp *interp,RECLASS reclass,int size);
  43 int DumpReclassToProgram(Tcl_Interp *interp,RASTER_OBJECT handle);
  44 /*
  45  * Cleans up global data, used by raster objects. Probably, should also
  46  * clean up all existing raster objects (?)
  47  */
  48
  49 EXPORT(void, Fgis_DeleteRaster)(ClientData table)
  50 {
  51 }
  52 /*
  53  * parses standard raster options -reclass and -table.
  54  * Returns pointer to reclass on success, NULL if option is not one of these
  55  * and FGIS_INVALID_RECLASS
  56  * if there was error parsing option argument.
  57  */
  58 RECLASS Fgis_StdRasterOptions(Tcl_Interp *interp, int size, char **arg)
  59 { RECLASS table;
  60    if (!strcmp(*arg,"-reclass")) {
  61       table=Fgis_ParseReclass(interp,*(arg+1),create_reclass_table(size),size);
  62    } else if (!strcmp(*arg,"-table")) {
  63       table=Fgis_ListToReclass(interp,*(arg+1),create_reclass_table(size),size);
  64    } else return NULL;
  65    if (!table)
  66      return FGIS_INVALID_RECLASS;
  67    else
  68      return table;
  69 }
  70
  71 /*
  72  * Implements raster Tcl command - various ways to create new raster object
  73  * (deals itself with existing raster files and calls CreateNewRaster for new
  74  * ones.
  75  */
  76
  77 EXPORT(int, Fgis_Raster)(ClientData data,Tcl_Interp *interp,
  78                    int argc,
  79                    char **argv)
  80 {
  81  RASTER_OBJECT curfile;
  82  XEPP* xptr;
  83  int reclass_size;
  84  EPP *(*openfunc)(char *)=open_epp;
  85     if (Fgis_CkArgs(interp,argc<2,argv[0],"?mode? filename ?options?")) {
  86         return TCL_ERROR;
  87     }
  88     if (!strcmp(argv[1],"load")) {
  89        argc--;
  90        argv++;
  91        openfunc=load_epp;
  92     } else if (!strcmp(argv[1],"new")){
  93        return Fgis_CreateNewRaster(data,interp,argc-2,argv+2);
  94     }
  95     if (Fgis_CkArgs(interp,argc!=2&&argc!=4,argv[0],"filename ?options?")) {
  96       return TCL_ERROR;
  97     }
  98     /* Open file, or get pointer to existing one */
  99     xptr=Fgis_OpenXEPP(argv[1],open_epp);
 100      if (!xptr) {
 101         if (map_error==ME_NO_FILE||map_error==ME_CREATE_ERROR) {
 102           Tcl_AppendResult(interp,"Cannot open ",argv[1],": ",
 103                 Tcl_PosixError(interp),NULL);
 104         } else {
 105           Tcl_AppendResult(interp,"Invalid epp file:",argv[1],NULL);
 106         }
 107         return TCL_ERROR;
 108      }
 109     /* Create reclass table with size of epp_max, if file is not loaded,
 110        or 1<<kind if it is loaded. 8-bit files always have
 111        256-element reclass table, becouse some old software can
 112        create epp files with non-correct max-min */
 113
 114      if (xptr->editable) {
 115        reclass_size=1<<(xptr->e->kind);
 116      } else {
 117        reclass_size=xptr->e->kind==8?256:epp_table_size(xptr->e);
 118      }
 119
 120     /* Allocate new raster object */
 121     curfile=malloc(sizeof(struct RASTER_OBJECT));
 122     curfile->file=xptr;
 123     if (argc>2) {
 124       /* If there is reclass option, parse it */
 125       curfile->reclass=Fgis_StdRasterOptions(interp,reclass_size,argv+2);
 126       if (!(curfile->reclass)||curfile->reclass==FGIS_INVALID_RECLASS) {
 127           curfile->reclass=NULL;
 128           Fgis_CloseXEPP(curfile->file);
 129           return TCL_ERROR;
 130       }
 131     } else {
 132       /* otherwise create empty reclass */
 133       curfile->reclass=create_reclass_table(reclass_size);
 134     }
 135      curfile->next=firstobject;
 136      firstobject=curfile;
 137      return Fgis_CreateObjectCommand(interp,"raster",Fgis_RasterObj,
 138             (ClientData)curfile, Fgis_DeleteRasterObj);
 139 }
 140
 141 /*
 142  * Creates new raster object. returns TCL_ERROR if file already exists
 143  *
 144  */
 145
 146 EXPORT(int, Fgis_CreateNewRaster)(ClientData data,Tcl_Interp *interp,int argc,
 147                          char **argv)
 148 { struct stat BUF;
 149   RECLASS table=NULL,t1;
 150   int i,width=0,height=0,offsite=-1,explicit_limits=0;
 151   EPP *e,*pattern=NULL;
 152   XEPP *xptr;
 153   double X1,Y1,X2,Y2;
 154   RASTER_OBJECT object;
 155   int bits=16;
 156   char filename[1024];
 157    strcpy(filename,argv[1]);
 158    /* checking file existence */
 159     if (stat(filename,&BUF)||errno!=ENOENT) {
 160         Tcl_AppendResult(interp,"File ",filename," already exists",NULL);
 161         return TCL_ERROR;
 162     }
 163  /* Option parsing loop */
 164     for(i=0;i<argc;i++) {
 165         if (!strcmp(argv[i],"-like")) {
 166             RASTER_OBJECT p;
 167             if (!(p=Fgis_GetRaster(interp,argv[++i])))
 168                 return TCL_ERROR;
 169             pattern=epp(p);
 170         } else if (!strcmp(argv[i],"-width")) {
 171             if (Fgis_GetInt(interp,argc,argv,++i,0,32767,&width,"raster width")
 172                    !=TCL_OK)
 173                 return TCL_ERROR;
 174         } else if (!strcmp(argv[i],"-height")) {
 175             if (Fgis_GetInt(interp,argc,argv,++i,0,32767,&height,
 176                     "raster height")!=TCL_OK)
 177                 return TCL_ERROR;
 178         } else if (!strcmp(argv[i],"-offsite")) {
 179             if (Fgis_GetInt(interp,argc,argv,++i,0,65535,&offsite,
 180                     "offsite value")!=TCL_OK)
 181                 return TCL_ERROR;
 182         } else if (!strcmp(argv[i],"-8bit")) {
 183             bits=8;
 184             if (table) {
 185               table=realloc(table,256*sizeof(short));
 186             }
 187         } else if (!strcmp(argv[i],"-limits")) {
 188             if (++i==argc)
 189                 ERROR_MESSAGE("List of four reals expected",TCL_STATIC);
 190             if (Fgis_GetLimits(interp,argv[i],&X1,&Y1,&X2,&Y2)==TCL_ERROR)
 191                 return TCL_ERROR;
 192             explicit_limits=1;
 193         } else if ((t1=Fgis_StdRasterOptions(interp,1<<bits,argv+i))) {
 194            if (table) {
 195              Tcl_SetResult(interp, "Duplicate reclass specification",
 196                     TCL_STATIC);
 197              free(table);
 198              if (t1!=FGIS_INVALID_RECLASS) free(t1);
 199              return TCL_ERROR;
 200            }
 201            if (t1==FGIS_INVALID_RECLASS) return TCL_ERROR;
 202            table=t1;
 203            i++;
 204         } else {
 205             if (table) free(table);
 206             ERROR_MESSAGE("Wrong option. Must be one of -like -width"
 207                   "-height -offsite -limits -8bit -reclass -table",TCL_STATIC);
 208         }
 209     }
 210      Create16bit=bits==16;
 211
 212     if (bits==8&&offsite>255&&offsite!=65535)  {
 213          if (table) free(table);
 214          ERROR_MESSAGE("Invalid offsite for 8-bit file",TCL_STATIC);
 215     }
 216
 217     if(offsite<0) {
 218       if (bits==8)
 219            offsite=255;
 220        else
 221          offsite=65535;
 222     }
 223     if (!pattern) {
 224         if (!width || !height) {
 225             if (table) free(table);
 226             ERROR_MESSAGE("Raster size is not specified",TCL_STATIC);
 227         }
 228         if (!explicit_limits||X1==X2||Y1==Y2) {
 229             X1=0;X2=width;Y1=height;Y2=0;
 230         }
 231         e=creat_epp(filename,1,1,width,height,X1,Y1,X2,Y2,100,0,offsite);
 232         if (!e) {
 233            if (table) free(table);
 234            Tcl_AppendResult(interp,"Error creating file ",filename,":",
 235              Tcl_PosixError(interp),NULL);
 236            return TCL_ERROR;
 237         }
 238     } else {
 239         if (width || height) {
 240              if (table) free(table);
 241             ERROR_MESSAGE("Raster size is specified twice",TCL_STATIC);
 242         }
 243         if (explicit_limits) {
 244             if (table) free(table);
 245             ERROR_MESSAGE("Coordinate system is specified twice",TCL_STATIC);
 246         }
 247         e=creat_epp_as(filename,pattern);
 248         if (!e) {
 249            if (table) free(table);
 250            Tcl_AppendResult(interp,"Error creating file ",filename,":",
 251              Tcl_PosixError(interp),NULL);
 252            return TCL_ERROR;
 253         }
 254         e->offsite=offsite;
 255     }
 256     /* Fill entire raster with offsite value */
 257     /* And load it into memory to make editable */
 258     load_new_epp(e);
 259     /* Now create raster object around it */
 260      xptr=Fgis_NewXEPP(e,filename);
 261      xptr->editable=1;
 262      object=malloc(sizeof (struct RASTER_OBJECT));
 263      object->file=xptr;
 264      if (!table) {
 265         object->reclass=create_reclass_table(1<<e->kind);
 266      } else {
 267         object->reclass=table;
 268      }
 269      object->next=firstobject;
 270      firstobject=object;
 271      return  Fgis_CreateObjectCommand(interp,"raster",Fgis_RasterObj,
 272           (ClientData)object, Fgis_DeleteRasterObj);
 273 }
 274
 275 /*
 276  * Fgis_GetRaster
 277  * returns pointer to raster object, associated with tcl command
 278  *
 279  *
 280  */
 281 EXPORT(RASTER_OBJECT, Fgis_GetRaster)(Tcl_Interp *interp, char *name)
 282 {
 283     Tcl_CmdInfo info;char buffer[255]="Invalid palette: ";
 284     if (!Tcl_GetCommandInfo(interp,name,&info)) return NULL;
 285     if (info.proc!=&Fgis_RasterObj) {
 286        strcat(buffer,name);
 287        Tcl_SetResult(interp,name,TCL_VOLATILE);
 288        return NULL;
 289     }
 290     return (RASTER_OBJECT)(info.clientData);
 291 }
 292
 293
 294 /*
 295  *
 296  *  Fgis_DeleteRasterObj - deletes raster object, associated with
 297  *  command. Called as Tcl_CmdDeleteProc, but not for base rasters.
 298  *
 299  */
 300
 301 EXPORT(void, Fgis_DeleteRasterObj)(ClientData data)
 302 { RASTER_OBJECT r=(RASTER_OBJECT)data;
 303   RASTER_OBJECT prev;
 304   /* Remove it from list */
 305   if (firstobject==r) {
 306       firstobject=r->next;
 307   } else {
 308     for (prev=firstobject;prev->next!=r;prev=prev->next);
 309     prev->next=r->next;
 310   }
 311   Fgis_CloseXEPP(r->file);
 312   free(r->reclass);
 313   free(r);
 314 }
 315
 316 /*
 317  * Fgis_RasterObj - implementation of object command for rasters.
 318  * mostly calls some procedure according to first argument. But few commands,
 319  * which returns info of object itself, rather then of raster or reclass,
 320  * are inlined, as well as delete command.
 321  */
 322 EXPORT(int, Fgis_RasterObj)(ClientData data,Tcl_Interp *interp, int argc, char **argv)
 323 {
 324     RASTER_OBJECT object=(RASTER_OBJECT)data;
 325
 326     if (Fgis_CkArgs(interp,argc<2,argv[0], " option ?args?"))
 327         return TCL_ERROR;
 328
 329     /* Checking for valid options*/
 330     /* delete */
 331     if (!strcmp(argv[1],"delete")) {
 332         if (Fgis_CkArgs(interp,argc!=2,argv[0], "delete"))
 333             return TCL_ERROR;
 334         Tcl_DeleteCommand(interp,argv[0]);
 335         return TCL_OK;
 336      /* filename */
 337     } else if (!strcmp(argv[1],"filename")) {
 338       if (Fgis_CkArgs(interp,argc!=2,argv[0],
 339             " filename"))
 340          return TCL_ERROR;
 341         RETURN(object->file->filename,TCL_VOLATILE);
 342      /*load */
 343     } else if (!strcmp(argv[1],"cell")) {
 344         return ReturnCell(object,interp,argc,argv);
 345     } else if (!strcmp(argv[1],"save")) {
 346         return SaveRaster(object,interp,argc,argv);
 347     } else if (!strcmp(argv[1],"bpp"))  {
 348         return ReturnBPP(object,interp,argc,argv);
 349     } else if (!strcmp(argv[1],"cache")) {
 350         return SetRasterCache(object,interp,argc,argv);
 351     }else if (!strcmp(argv[1],"get")) {
 352         return ReturnCellValue(object,interp,argc,argv);
 353     }else if (!strcmp(argv[1],"limits")) {
 354         return ReturnLimits(object,interp,argc,argv);
 355     }else if (!strcmp(argv[1],"comment")) {
 356         return RasterComments(object,interp,argc,argv);
 357     }else if (!strcmp(argv[1],"reclass")) {
 358         return ChangeReclass(object,interp,argc,argv);
 359     }else if (!strcmp(argv[1],"offsite")) {
 360         return OffsiteValue(object,interp,argc,argv);
 361     }else if (!strcmp(argv[1],"row")) {
 362         return ReturnRow(object,interp,argc,argv);
 363     }else if (!strcmp(argv[1],"col")) {
 364         return ReturnCol(object,interp,argc,argv);
 365     }else if (!strcmp(argv[1],"put")) {
 366         return ModifyCell(object,interp,argc,argv);
 367     }else if (!strcmp(argv[1],"line")) {
 368         return PlotLine(object,interp,argc,argv);
 369     }else if (!strcmp(argv[1],"frame")) {
 370         return PlotFrame(object,interp,argc,argv);
 371     }else if (!strcmp(argv[1],"box")) {
 372         return PlotBox(object,interp,argc,argv);
 373     }else if (!strcmp(argv[1],"fill")) {
 374         return FillContour(object,interp,argc,argv);
 375     }else if (!strcmp(argv[1],"circle")) {
 376         return PlotCircle(object,interp,argc,argv);
 377     }else if (!strcmp(argv[1],"classes")) {
 378         return ListClasses(object,interp,argc,argv);
 379     }else if (!strcmp(argv[1],"count")) {
 380         return CountClass(object,interp,argc,argv);
 381     }else if (!strcmp(argv[1],"unused")) {
 382         return FindUnused(object,interp,argc,argv);
 383     }else if (!strcmp(argv[1],"extents")) {
 384         return CalcExtents(object,interp,argc,argv);
 385     }else if (!strcmp(argv[1],"celllimits")) {
 386         return ReturnCellLimits(object,interp,argc,argv);
 387     }else if (!strcmp(argv[1],"xleft")) {
 388         return XLeft(object,interp,argc,argv);
 389     }else if (!strcmp(argv[1],"xright")) {
 390         return XRight(object,interp,argc,argv);
 391     }else if (!strcmp(argv[1],"ytop")) {
 392         return YTop(object,interp,argc,argv);
 393     }else if (!strcmp(argv[1],"ybottom")) {
 394         return YBottom(object,interp,argc,argv);
 395     }else if (!strcmp(argv[1],"shift")) {
 396         return ShiftCoords(object,interp,argc,argv);
 397     }else if (!strcmp(argv[1],"max")) {
 398         return GetMaxValue(object,interp,argc,argv);
 399     }else if (!strcmp(argv[1],"min")) {
 400         return GetMinValue(object,interp,argc,argv);
 401     }else if (!strcmp(argv[1],"unit")) {
 402         return CellUnit(object,interp,argc,argv);
 403     }else {
 404         Tcl_SetResult(interp,"Invalid option. Should be one of "
 405                     "box bpp cache circle classes col comment count delete "
 406                     "extents filename fill first_row first_col frame get "
 407                     "limits line load max min "
 408                     "polygon put reclass row save shift transect unused unit"
 409                     "xleft xright ybottom ytop",TCL_STATIC);
 410         return TCL_ERROR;
 411     }
 412 }
 413 int SaveRaster (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 414   NOT_YET
 415
 416
 417 int ReturnCellValue (RASTER_OBJECT data, Tcl_Interp *interp, int argc,
 418    char **argv)
 419 {   char buffer[10];
 420     EPP *e=data->file->e;
 421     double X,Y;
 422     int value;
 423
 424     if (Fgis_CkArgs(interp,(argc != 4&&!(argc==5&&!strcmp(argv[4],"-base"))),
 425          argv[0],"get x y ?-base?")) return TCL_ERROR;
 426     if (Tcl_GetDouble(interp,argv[2],&X)!=TCL_OK)
 427         return TCL_ERROR;
 428     if (Tcl_GetDouble(interp,argv[3],&Y)!=TCL_OK)
 429         return TCL_ERROR;
 430     value= epp_get(e,epp_col(e,X),epp_row(e,Y));
 431     sprintf(buffer,"%d",argc!=5?data->reclass[value]:value);
 432     RETURN(buffer,TCL_VOLATILE);
 433 }
 434
 435 int ReturnLimits (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 436 {
 437     char buffer[255];
 438     EPP *e=epp(data);
 439
 440     if (Fgis_CkArgs(interp,argc!=2,argv[0],"limits"))
 441         return TCL_ERROR;
 442     Tcl_PrintDouble(interp,e->XLeft,buffer);
 443     Tcl_AppendElement(interp,buffer);
 444     Tcl_PrintDouble(interp,e->YBottom,buffer);
 445     Tcl_AppendElement(interp,buffer);
 446     Tcl_PrintDouble(interp,e->XRight,buffer);
 447     Tcl_AppendElement(interp,buffer);
 448     Tcl_PrintDouble(interp,e->YTop,buffer);
 449     Tcl_AppendElement(interp,buffer);
 450     return TCL_OK;
 451 }
 452 int RasterComments (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 453 {
 454  EPP *e=epp(data);
 455
 456  switch (argc) {
 457  case 3: /* set comment */
 458         if (!data->file->editable)
 459           ERROR_MESSAGE("Raster is read-only.",TCL_STATIC);
 460         setcomment(e,argv[2]);
 461          /* NO BREAK HERE */
 462  case 2: /* return comment */
 463         RETURN(getcomment(e),TCL_VOLATILE);
 464  default:{
 465           Tcl_AppendResult(interp,"Wrong #args. Should be ",
 466                  argv[0], " comment ?string?",NULL);
 467           return TCL_ERROR;
 468          }
 469  }
 470 }
 471 int OffsiteValue (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 472 {
 473     EPP *e=epp(data);
 474     char result[7];
 475     int new_offsite;
 476
 477     switch (argc) {
 478     case 3: /*set new offsite value*/
 479                if (!data->file->editable);
 480                        ERROR_MESSAGE("Raster is read-only",TCL_STATIC);
 481                if (Fgis_GetInt(interp,argc,argv,2,-1,
 482                        (1<<e->kind)-1,&new_offsite,"offsite value")!=TCL_OK)
 483                    return TCL_ERROR;
 484                if (new_offsite<0) new_offsite=65535;
 485                e->offsite=new_offsite;
 486               /* NO BREAK HERE!!*/
 487    case 2: /* return offsite value */
 488            sprintf(result,"%d",data->reclass[e->offsite]);
 489            RETURN(result,TCL_VOLATILE);
 490   default: Tcl_AppendResult(interp,
 491             "Wrong # args. Should be ",argv[0]," offsite ?value?",NULL);
 492            return TCL_ERROR;
 493   }
 494 }
 495
 496 int CellUnit (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 497 {
 498     EPP *e=epp(data);
 499     static char *units[]={"undefined","ft","m","km","mile","ha","acre",NULL};
 500     EPPHEADER h;
 501
 502     get_epp_header(e,&h);
 503     if (h.area_unit>6) h.area_unit=0;
 504     if (argc==3) {
 505        char **u;
 506        if (!data->file->editable)
 507            ERROR_MESSAGE("Raster is read-only",TCL_STATIC);
 508        for(u=units;*u&&strcmp(*u,argv[2]);u++);
 509        if (!*u) {
 510            Tcl_SetResult(interp,"Wrong area unit. should be one of: ",
 511                      TCL_STATIC);
 512            for(u=units;*u;u++)
 513                  Tcl_AppendResult(interp,*u," ",NULL);
 514            return TCL_ERROR;
 515        }
 516        h.area_unit=u-units;
 517        change_epp_header(e,h);
 518     } else if (Fgis_CkArgs(interp,argc!=2,argv[0],
 519                 " unit ?new-unit?")) {
 520         return TCL_ERROR;
 521     }
 522     RETURN(units[(int)h.area_unit],TCL_STATIC);
 523 }
 524
 525 int GetMaxValue (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 526 {
 527     char result[255];
 528
 529     if (Fgis_CkArgs(interp,argc!=2,argv[0],"max"))
 530        return TCL_ERROR;
 531     sprintf(result,"%d",Fgis_RasterMax(data));
 532     RETURN(result,TCL_VOLATILE);
 533 }
 534
 535 int GetMinValue (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 536 {
 537     char result[255];
 538
 539     if (Fgis_CkArgs(interp,argc!=2,argv[0],"min"))
 540        return TCL_ERROR;
 541     sprintf(result,"%d",Fgis_RasterMin(data));
 542     RETURN(result,TCL_VOLATILE);
 543 }
 544
 545 int ReturnBPP (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 546 {
 547     EPP *e=epp(data);
 548     if (Fgis_CkArgs(interp,argc!=2,argv[0],"bpp"))
 549        return TCL_ERROR;
 550     RETURN(e->kind==8?"8":"16",TCL_STATIC);
 551 }
 552
 553 int ReturnCell (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 554 {
 555   EPP *e=epp(data);
 556   char result[32];
 557   if (Fgis_CkArgs(interp,argc!=2&&!(argc==3&&!strcmp(argv[2],"-area")),argv[0],
 558           "cell ?-area?"))
 559        return TCL_ERROR;
 560   if (argc==3) {
 561      Tcl_PrintDouble(interp,e->cell_area,result);
 562   } else {
 563      Tcl_PrintDouble(interp,fabs(e->XRight-e->XLeft)/(e->lc-e->fc),result);
 564   }
 565   RETURN(result,TCL_VOLATILE);
 566 }
 567 int ReturnCellLimits (RASTER_OBJECT data, Tcl_Interp *interp,
 568        int argc, char **argv)
 569 {
 570   EPP *e=epp(data);
 571   char result[7];
 572   if (Fgis_CkArgs(interp,argc!=2,argv[0],"celllimits"))
 573     return TCL_ERROR;
 574   sprintf(result,"%d",e->fc);
 575   Tcl_AppendElement(interp,result);
 576   sprintf(result,"%d",e->lr-1);
 577   Tcl_AppendElement(interp,result);
 578   sprintf(result,"%d",e->lc-1);
 579   Tcl_AppendElement(interp,result);
 580   sprintf(result,"%d",e->fr);
 581   Tcl_AppendElement(interp,result);
 582   return TCL_OK;
 583 }
 584
 585 int XLeft (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 586 {
 587    EPP *e=epp(data);
 588    char result[64];
 589    if (argc==3) {
 590       EPPHEADER h;
 591       double tmp;
 592
 593       if (!data->file->editable)
 594           ERROR_MESSAGE("Raster is read-only",TCL_STATIC);
 595       if (Tcl_GetDouble(interp,argv[2],&tmp)!=TCL_OK)
 596           return TCL_ERROR;
 597       get_epp_header(e,&h);
 598       h.fcx=tmp;
 599       e->XLeft=tmp;
 600       change_epp_header(e,h);
 601    } else if (argc!=2) {
 602      Tcl_AppendResult(interp,"Wrong # args. Should be", argv[0]," ",
 603              argv[1]," ?value?", NULL);
 604  }
 605  Tcl_PrintDouble(interp,e->XLeft,result);
 606  RETURN(result,TCL_VOLATILE);
 607 }
 608 int XRight (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 609 {
 610    EPP *e=epp(data);
 611    char result[64];
 612    if (argc==3) {
 613       EPPHEADER h;
 614       double tmp;
 615
 616       if (!data->file->editable)
 617           ERROR_MESSAGE("Raster is read-only",TCL_STATIC);
 618       if (Tcl_GetDouble(interp,argv[2],&tmp)!=TCL_OK)
 619           return TCL_ERROR;
 620       get_epp_header(e,&h);
 621       h.lcx=tmp;
 622       e->XRight=tmp;
 623       change_epp_header(e,h);
 624    } else if (argc!=2) {
 625      Tcl_AppendResult(interp,"Wrong # args. Should be", argv[0]," ",
 626              argv[1]," ?value?", NULL);
 627  }
 628  Tcl_PrintDouble(interp,e->XRight,result);
 629  RETURN(result,TCL_VOLATILE);
 630 }
 631 int YTop (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 632 {
 633    EPP *e=epp(data);
 634    char result[64];
 635    if (argc==3) {
 636       EPPHEADER h;
 637       double tmp;
 638
 639       if (!data->file->editable)
 640           ERROR_MESSAGE("Raster is read-only",TCL_STATIC);
 641       if (Tcl_GetDouble(interp,argv[2],&tmp)!=TCL_OK)
 642           return TCL_ERROR;
 643       get_epp_header(e,&h);
 644       h.fry=tmp;
 645       e->YTop=tmp;
 646       change_epp_header(e,h);
 647    } else if (argc!=2) {
 648      Tcl_AppendResult(interp,"Wrong # args. Should be", argv[0]," ",
 649              argv[1]," ?value?", NULL);
 650  }
 651  Tcl_PrintDouble(interp,e->YTop,result);
 652  RETURN(result,TCL_VOLATILE);
 653 }
 654 int YBottom (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 655 {
 656    EPP *e=epp(data);
 657    char result[64];
 658    if (argc==3) {
 659       EPPHEADER h;
 660       double tmp;
 661
 662       if (!data->file->editable)
 663           ERROR_MESSAGE("Raster is read-only",TCL_STATIC);
 664       if (Tcl_GetDouble(interp,argv[2],&tmp)!=TCL_OK)
 665           return TCL_ERROR;
 666       get_epp_header(e,&h);
 667       h.lry=tmp;
 668       e->YBottom=tmp;
 669       change_epp_header(e,h);
 670    } else if (argc!=2) {
 671      Tcl_AppendResult(interp,"Wrong # args. Should be", argv[0]," ",
 672              argv[1]," ?value?", NULL);
 673  }
 674  Tcl_PrintDouble(interp,e->YBottom,result);
 675  RETURN(result,TCL_VOLATILE);
 676 }
 677
 678 int ShiftCoords (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 679 {
 680     EPP *e=epp(data);
 681     EPPHEADER h;
 682     double dx,dy;
 683     if (Fgis_CkArgs(interp,argc!=4,argv[0], "shift dx dy")) {
 684       return TCL_ERROR;
 685     }
 686
 687     if (!data->file->editable)
 688         ERROR_MESSAGE("Raster is read-only",TCL_STATIC);
 689     if (Tcl_GetDouble(interp,argv[2],&dx)!=TCL_OK)
 690        return TCL_ERROR;
 691     if (Tcl_GetDouble(interp,argv[3],&dy)!=TCL_OK)
 692        return TCL_ERROR;
 693     get_epp_header(e,&h);
 694     h.fcx=e->XLeft+=dx;
 695     h.lcx=e->XRight+=dx;
 696     h.fry=e->YTop+=dy;
 697     h.lry=e->YBottom+=dy;
 698     change_epp_header(e,h);
 699     return TCL_OK;
 700 }
 701
 702 int ChangeReclass(RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 703 {
 704     RECLASS tmp;
 705     EPP *e=epp(data);
 706     if (argc==2)
 707       return DumpReclassToList(interp,data->reclass,epp_table_size(e));
 708     if (!strcmp(argv[2],"-table")) {
 709         if (Fgis_CkArgs(interp,(argc!=4), argv[0],"reclass -table list")) {
 710            return TCL_ERROR;
 711         }
 712         tmp=Fgis_ListToReclass(interp,argv[3],data->reclass,
 713                 epp_table_size(e));
 714     } else
 715     if (!strcmp(argv[2],"-statements")) {
 716         if (Fgis_CkArgs(interp,argc!=3, argv[0], "reclass -statements")) {
 717             return TCL_ERROR;
 718         }
 719         return DumpReclassToProgram(interp,data);
 720     } else {
 721         if (Fgis_CkArgs(interp,argc!=3, argv[0], " reclass ?option? arg")) {
 722             return TCL_ERROR;
 723         }
 724         tmp=Fgis_ParseReclass(interp,argv[2],data->reclass,epp_table_size(e));
 725     }
 726     if(!tmp) return TCL_ERROR;
 727     free(data->reclass);
 728     data->reclass=tmp;
 729     return TCL_OK;
 730 }
 731
 732 int ReturnRow (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 733 {
 734     EPP *e=epp(data);
 735     char result[7];
 736     int row;
 737     double Y;
 738     if (Fgis_CkArgs(interp,argc != 3, argv[0],"row y")) {
 739         return TCL_ERROR;
 740     }
 741     if (Tcl_GetDouble(interp,argv[2],&Y)!=TCL_OK) return TCL_ERROR;
 742     row=epp_row(e,Y);
 743     if (row<e->fr||row>=e->lr)
 744        ERROR_MESSAGE("Y coordinate outside file boundaries",TCL_STATIC);
 745     sprintf(result,"%d",row);
 746     RETURN(result,TCL_VOLATILE);
 747 }
 748 int ReturnCol (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 749 {
 750     EPP *e=epp(data);
 751     char result[7];
 752     int col;
 753     double X;
 754     if (Fgis_CkArgs(interp,argc != 3,argv[0],"col x")) {
 755         return TCL_ERROR;
 756     }
 757     if (Tcl_GetDouble(interp,argv[2],&X)!=TCL_OK) return TCL_ERROR;
 758     col=epp_col(e,X);
 759     if (col<e->fc||col>=e->lc)
 760        ERROR_MESSAGE("X coordinate outside file boundaries",TCL_STATIC);
 761     sprintf(result,"%d",col);
 762     RETURN(result,TCL_VOLATILE);
 763 }
 764
 765 int SetRasterCache (RASTER_OBJECT data, Tcl_Interp *interp, int argc,
 766     char **argv)
 767 {
 768 EPP *e=epp(data);
 769 char result[7];
 770
 771 if (Fgis_CkArgs(interp,argc>3,argv[0],"cache ?cache-size?")) {
 772    return TCL_ERROR;
 773 }
 774 if (argc==2) {
 775     if (e->mode&MAP_LOADED)
 776       RETURN("loaded",TCL_STATIC);
 777    /* Otherwise cache size would be returned, as it is done after change */
 778 } else {
 779    int newsize;
 780    if (e->mode&MAP_LOADED) {
 781       Tcl_AppendResult(interp,argv[0]," already loaded into memory",NULL);
 782       return TCL_ERROR;
 783    }
 784    if (Fgis_GetInt(interp,argc,argv,2,-1,CACHE_THRESHOLD/
 785            (e->width*(e->kind>>3)),
 786            &newsize,"cache size")!=TCL_OK) return TCL_ERROR;
 787    set_epp_cache(e,newsize);
 788 }
 789 sprintf(result,"%d",e->cache_size);
 790 RETURN(result,TCL_VOLATILE);
 791 }
 792 /*
 793 # $raster put value x y
 794 # modifies raster (if modifiable). Value should be specified in
 795 # terms of base raster, not of current raster object.
 796 #
 797 */
 798 int ModifyCell (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 799 { EPP *e=epp(data);
 800   int row,col,value;
 801   double X,Y;
 802   if (Fgis_CkArgs(interp,argc!=5,argv[0]," put value x y")) {
 803     return TCL_ERROR;
 804   }
 805   if (!(e->mode&MAP_LOADED)) {
 806      Tcl_AppendResult(interp,"Raster \"",argv[0],"\" is read-only",NULL);
 807      return TCL_ERROR;
 808   }
 809   if (Tcl_GetInt(interp,argv[2],&value)==TCL_ERROR) {
 810      return TCL_ERROR;
 811   }
 812   if (Tcl_GetDouble(interp,argv[3],&X)==TCL_ERROR) {
 813      return TCL_ERROR;
 814   }
 815   if (Tcl_GetDouble(interp,argv[4],&Y)==TCL_ERROR) {
 816      return TCL_ERROR;
 817   }
 818   if (value<0||value>=MAX_EPP_CLASS) {
 819     Tcl_AppendResult(interp,"Value ",argv[2]," is out of range",NULL);
 820     return TCL_ERROR;
 821   }
 822   row=epp_row(e,Y);
 823   col=epp_col(e,X);
 824   if (row<e->fr||row>=e->lr||col<e->fc||col>=e->lc) {
 825     Tcl_AppendResult(interp,"Point ",argv[3],",",argv[4]," is out of physical"
 826            " boundaries of file", NULL);
 827     return TCL_ERROR;
 828   }
 829   /* If we have value, which does't fit to 8 */
 830   if (value>=(1<<e->kind)) {
 831      Tcl_AppendResult(interp,"Value ",argv[2]," doesn't fit into raster ",
 832            argv[0],NULL);
 833      return TCL_ERROR;
 834   }
 835   epp_put(e,col,row,value);
 836   return TCL_OK;
 837 }
 838 int PlotLine (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 839     NOT_YET
 840 int PlotFrame (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 841     NOT_YET
 842 int PlotBox (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 843     NOT_YET
 844 int FillContour (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 845    NOT_YET
 846 int PlotCircle (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 847    NOT_YET
 848
 849 int ListClasses (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 850 {
 851     int value,i;
 852     char result[7];
 853     EPP *e=epp(data);
 854     if( Fgis_CkArgs(interp,argc!=3,argv[0]," classes value")) {
 855        return TCL_ERROR;
 856     }
 857     if (Fgis_GetInt(interp,argc,argv,2,0,65535,&value,"cell value")!=TCL_OK)
 858         return TCL_ERROR;
 859
 860     for(i=e->min;i<=e->max;i++)
 861         if ((i!=e->offsite)&&(data->reclass[i]==value)) {
 862             sprintf(result,"%d",i);
 863             Tcl_AppendElement(interp,result);
 864         }
 865     return TCL_OK;
 866 }
 867 int CountClass (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 868   NOT_YET
 869 int FindUnused (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 870   NOT_YET
 871 int CalcExtents (RASTER_OBJECT data, Tcl_Interp *interp, int argc, char **argv)
 872    NOT_YET
 873
 874 #define LIST_ERROR(msg,mode) { Tcl_Free((char *)listv);free(table);\
 875                 Tcl_SetResult(interp,msg,mode);\
 876                 return NULL;}
 877
 878 RECLASS Fgis_ListToReclass(Tcl_Interp *interp,char *list,RECLASS src,int size)
 879 {
 880     int listc,i,index,value;
 881     char **listv;
 882     RECLASS table;
 883
 884     if (Tcl_SplitList(interp,list,&listc,&listv)!=TCL_OK)
 885         return NULL;
 886     table=memcpy(malloc((size+1)*sizeof(short int)),
 887             src,
 888             (size+1)*sizeof(short int));
 889
 890     for(i=0;i<listc;i++) {
 891         int pairc; char **pairv;
 892         if (Tcl_SplitList(interp,listv[i],&pairc,&pairv)!=TCL_OK||pairc!=2)
 893             LIST_ERROR("Each element of list shold be pair of integers",
 894                     TCL_STATIC);
 895         if (Tcl_GetInt(interp,pairv[0],&index)!=TCL_OK||
 896                 Tcl_GetInt(interp,pairv[1],&value)!=TCL_OK)
 897             LIST_ERROR("Integer value expected",TCL_STATIC);
 898         if (index<=size&&index>=0)
 899             table[index]=value;
 900     }
 901     return table;
 902 }
 903
 904 unsigned char *reclass_program;
 905
 906 static int get_prog_char()
 907 {
 908     return *(reclass_program++);
 909 }
 910
 911 RECLASS Fgis_ParseReclass(Tcl_Interp *interp,char *program,RECLASS src,int size)
 912 {
 913     RECLASS table;
 914     reclass_program = (unsigned char *) program;
 915     table=memcpy(malloc((size+1)*sizeof(short int)),
 916             src,
 917             (size+1)*sizeof(short int));
 918     table = parse_statements(size,table, get_prog_char);
 919     if (!table) {
 920         Tcl_SetResult (interp, "Error in reclass statements", TCL_STATIC);
 921     }
 922     return table;
 923 }
 924
 925
 926
 927 int DumpReclassToList(Tcl_Interp *interp,RECLASS reclass,int size)
 928 {
 929     int i;
 930     char result[15];
 931
 932     for (i=0;i<=size;i++) {
 933         if (reclass[i]!=i) {
 934             sprintf(result,"%d %d",i,reclass[i]);
 935             Tcl_AppendElement(interp,result);
 936         }
 937     }
 938     return TCL_OK;
 939 }
 940 struct STATEMENT {int newclass;
 941                   char *list;
 942                   int lastclass;
 943                   int rangestart;
 944                  };
 945
 946 char *addtostring(char *source,char separator,int class)
 947 {
 948     char tmp[10];
 949     sprintf(tmp,"%c%d",separator,class);
 950     if(!source)
 951         return strcpy(malloc(strlen(tmp)+1),tmp);
 952     else
 953         return strcat(realloc(source,strlen(source)+strlen(tmp)+1),tmp);
 954
 955 }
 956
 957 void add_class(struct STATEMENT *stmt,int class)
 958 {
 959     if (!(stmt->list)) {
 960         stmt->list=addtostring(NULL,'=',class);
 961         stmt->rangestart=stmt->lastclass=class;
 962     } else
 963     if ((stmt->lastclass==class-1)||
 964             ((class==stmt->newclass+1)&&(stmt->lastclass==class-2))) {
 965         stmt->lastclass=class;
 966     } else {
 967         if (stmt->lastclass!=stmt->rangestart)
 968             stmt->list=addtostring(stmt->list,':',stmt->lastclass);
 969         stmt->list=addtostring(stmt->list,' ',class);
 970         stmt->rangestart=stmt->lastclass=class;
 971     }
 972 }
 973 int DumpReclassToProgram(Tcl_Interp *interp,RASTER_OBJECT data)
 974 {
 975     int max=Fgis_RasterMax(data)+1,i;
 976     struct STATEMENT *p,*cur;
 977     EPP *e=epp(data);
 978     RECLASS r=data->reclass;
 979     p=malloc(max*sizeof(struct STATEMENT));
 980     memset(p,0,max*sizeof(struct STATEMENT));
 981     for(i=0;i<max;i++)p[i].newclass=i;
 982     for(i=0;i<=e->max;i++) {
 983         if(i!=e->offsite&&r[i]!=i)
 984             add_class(p+r[i],i);
 985     }
 986     for(i=0,cur=p;i<max;i++,cur++) {
 987         if (cur->list) {
 988             char n[7];
 989             if (cur->lastclass!=cur->rangestart)
 990                 cur->list=addtostring(cur->list,':',cur->lastclass);
 991             sprintf(n,"%d",cur->newclass);
 992             Tcl_AppendResult(interp,n,cur->list,"\n",NULL);
 993             free(cur->list);
 994        }
 995     }
 996     free(p);
 997     return TCL_OK;
 998 }