shputils.c

Source code, and some other odds and ends can be downloaded from http://shapelib.maptools.org/dl.

		mean=dsum/maxrec;
		sprintf(stmp,"%%.%df to %%.%df \t(%%.%df)",iDecimals,iDecimals,iDecimals);
		if (dlow < dhigh)       printf(stmp,dlow,dhigh,mean);
		else if (dlow == dhigh) {
				        sprintf(stmp,"= %%.%df",iDecimals);
		                        printf(stmp,dlow);
		                        }
		else printf("No Values");
		break;

	    }

        }
	printf("\n");
}

int selectrec()
{
long int value, ty;

   ty = DBFGetFieldInfo( hDBF, iselectitem, NULL, &iWidth, &iDecimals);
      switch(ty)
      {
      case FTString:
        puts("Invalid Item");
        iselect=FALSE;
	break;
      case FTInteger:
        value = DBFReadIntegerAttribute( hDBF, iRecord, iselectitem );
        for (j = 0; j<selcount; j++)
          {
          if (selectvalues[j] == value)
               if (iunselect) return(0);  /* Keep this record */
                        else  return(1);  /* Skip this record */
          }
	break;
      case FTDouble:
        puts("Invalid Item");
        iselect=FALSE;
        break;
      }
      if (iunselect) return(1);  /* Skip this record */
               else  return(0);  /* Keep this record */
}


int check_theme_bnd()
{
    if ( (adfBoundsMin[0] >= cxmin) && (adfBoundsMax[0] <= cxmax) &&
         (adfBoundsMin[1] >= cymin) && (adfBoundsMax[1] <= cymax) )
    {   /** Theme is totally inside clip area **/
        if (ierase) nEntities=0; /** SKIP THEME  **/
        else   iclip=FALSE; /** WRITE THEME (Clip not needed) **/
    }
            
    if ( ( (adfBoundsMin[0] < cxmin) && (adfBoundsMax[0] < cxmin) ) ||
         ( (adfBoundsMin[1] < cymin) && (adfBoundsMax[1] < cymin) ) ||
         ( (adfBoundsMin[0] > cxmax) && (adfBoundsMax[0] > cxmax) ) ||
         ( (adfBoundsMin[1] > cymax) && (adfBoundsMax[1] > cymax) ) )
    {   /** Theme is totally outside clip area **/
        if (ierase) iclip=FALSE; /** WRITE THEME (Clip not needed) **/
             else   nEntities=0; /** SKIP THEME  **/
    }
            
    if (nEntities == 0)
        puts("WARNING: Theme is outside the clip area."); /** SKIP THEME  **/
}

clip_boundary()
{
    int  inside;
    int  prev_outside;
    int  i2;
    int  j2;
    
       /*** FIRST check the boundary of the feature ***/
       if ( ( (psCShape->dfXMin < cxmin) && (psCShape->dfXMax < cxmin) ) ||
            ( (psCShape->dfYMin < cymin) && (psCShape->dfYMax < cymin) ) ||
            ( (psCShape->dfXMin > cxmax) && (psCShape->dfXMax > cxmax) ) ||
            ( (psCShape->dfYMin > cymax) && (psCShape->dfYMax > cymax) ) )
            {   /** Feature is totally outside clip area **/
            	if (ierase) return(1); /** WRITE RECORD **/
            	     else   return(0); /** SKIP  RECORD **/
            }
       
       if ( (psCShape->dfXMin >= cxmin) && (psCShape->dfXMax <= cxmax) &&
            (psCShape->dfYMin >= cymin) && (psCShape->dfYMax <= cymax) )
            {   /** Feature is totally inside clip area **/
            	if (ierase) return(0); /** SKIP  RECORD **/
            	     else   return(1); /** WRITE RECORD **/
            }
            
       if (iinside) 
            { /** INSIDE * Feature might touch the boundary or could be outside **/
            if (ierase)  return(1); /** WRITE RECORD **/
                 else    return(0); /** SKIP  RECORD **/
            }
       
       if (itouch)
          {   /** TOUCH **/
          if ( ( (psCShape->dfXMin <= cxmin) || (psCShape->dfXMax >= cxmax) ) && 
                 (psCShape->dfYMin >= cymin) && (psCShape->dfYMax <= cymax)    )
               {   /** Feature intersects the clip boundary only on the X axis **/
               if (ierase) return(0); /** SKIP  RECORD **/
                    else   return(1); /** WRITE RECORD **/
               }

          if (   (psCShape->dfXMin >= cxmin) && (psCShape->dfXMax <= cxmax)   && 
               ( (psCShape->dfYMin <= cymin) || (psCShape->dfYMax >= cymax) )  )
               {   /** Feature intersects the clip boundary only on the Y axis **/
               if (ierase) return(0); /** SKIP  RECORD **/
                    else   return(1); /** WRITE RECORD **/
               }
               
          for( j2 = 0; j2 < psCShape->nVertices; j2++ ) 
               {   /** At least one vertex must be inside the clip boundary **/
               if ( (psCShape->padfX[j2] >= cxmin  &&  psCShape->padfX[j2] <= cxmax) ||
                    (psCShape->padfY[j2] >= cymin  &&  psCShape->padfY[j2] <= cymax)  )
                    if (ierase) return(0); /** SKIP  RECORD **/
                         else   return(1); /** WRITE RECORD **/
               }
               
          /** All vertices are outside the clip boundary **/ 
          if (ierase) return(1); /** WRITE RECORD **/
               else   return(0); /** SKIP  RECORD **/
          }   /** End TOUCH **/
          
       if (icut)
          {   /** CUT **/
          /*** Check each vertex in the feature with the Boundary and "CUT" ***/
          /*** THIS CODE WAS NOT COMPLETED!  READ NOTE AT THE BOTTOM ***/
          i2=0;
          prev_outside=FALSE;
          for( j2 = 0; j2 < psCShape->nVertices; j2++ ) 
             {
             inside = psCShape->padfX[j2] >= cxmin  &&  psCShape->padfX[j2] <= cxmax  &&
                      psCShape->padfY[j2] >= cymin  &&  psCShape->padfY[j2] <= cymax ;
                      
             if (ierase) inside=(! inside);
             if (inside)
                 {
                 if (i2 != j2)
                     {
                     if (prev_outside)
                         {
                         /*** AddIntersection(i2);   /*** Add intersection ***/
                         prev_outside=FALSE;
                         }
                     psCShape->padfX[i2]=psCShape->padfX[j2];     /** move vertex **/
                     psCShape->padfY[i2]=psCShape->padfY[j2];
                     }
                 i2++;
                 } else {
                 if ( (! prev_outside) && (j2 > 0) )
                     {
                     /*** AddIntersection(i2);   /*** Add intersection (Watch out for j2==i2-1) ***/
                     /*** Also a polygon may overlap twice and will split into a several parts  ***/
                     prev_outside=TRUE;
                     }
                 }
             }
             
        printf("Vertices:%d   OUT:%d   Number of Parts:%d\n",
                psCShape->nVertices,i2, psCShape->nParts );
               
             psCShape->nVertices = i2;
             
             if (i2 < 2) return(0); /** SKIP RECORD **/
             /*** (WE ARE NOT CREATING INTERESECTIONS and some lines could be reduced to one point) **/
        
             if (i2 == 0) return(0); /** SKIP  RECORD **/
                  else    return(1); /** WRITE RECORD **/
          }  /** End CUT **/
}


/************************************************************************/
/*                            strncasecmp2()                            */
/*                                                                      */
/*      Compare two strings up to n characters                          */
/*      If n=0 then s1 and s2 must be an exact match                    */
/************************************************************************/

int strncasecmp2(char *s1, char *s2, int n)

{
int j,i;
   if (n<1) n=strlen(s1)+1;
   for (i=0; i<n; i++)
   {
      if (*s1 != *s2)
      {
         if (*s1 >= 'a' && *s1 <= 'z') {
            j=*s1-32;
            if (j != *s2) return(*s1-*s2);
         } else {
            if (*s1 >= 'A' && *s1 <= 'Z') { j=*s1+32; }
                                   else   { j=*s1;    }
            if (j != *s2) return(*s1-*s2); 
         }
      }
      s1++;
      s2++;
   }
   return(0);
}


#define  NKEYS (sizeof(unitkeytab) / sizeof(struct unitkey))
findunit(unit)
   char *unit;
   {
   struct unitkey {
     char   *name;
     double value;
   } unitkeytab[] = {
     "CM",            39.37,
     "CENTIMETER",    39.37,
     "CENTIMETERS",   39.37,  /** # of inches * 100 in unit **/
     "METER",          3937,
     "METERS",         3937,
     "KM",          3937000,
     "KILOMETER",   3937000, 
     "KILOMETERS",  3937000,
     "INCH",            100,
     "INCHES",          100,
     "FEET",           1200,
     "FOOT",           1200,
     "YARD",           3600,
     "YARDS",          3600,       
     "MILE",        6336000,
     "MILES",       6336000  
   };

   double unitfactor=0;
   for (j = 0; j < NKEYS; j++) {
    if (strncasecmp2(unit, unitkeytab[j].name, 0) == 0) unitfactor=unitkeytab[j].value;
   }
   return(unitfactor);
}

/* -------------------------------------------------------------------- */
/*      Display a usage message.                                        */
/* -------------------------------------------------------------------- */
void error()
    {	
	puts( "The program will append to an existing shape file or it will" );
	puts( "create a new file if needed." );
	puts( "Only the items in the first output file will be preserved." );
	puts( "When an item does not match with the append theme then the item");
	puts( "might be placed to an existing item at the same position and type." );
	puts( "  OTHER FUNCTIONS:" );
	puts( "  - Describe all items in the dbase file (Use ALL for more than 5000 recs.)");
	puts( "  - Select a group of shapes from a comma separated selection list.");
	puts( "  - UnSelect a group of shapes from a comma separated selection list.");
	puts( "  - Clip boundary extent or by theme boundary." );
	puts( "      Touch writes all the shapes that touch the boundary.");
	puts( "      Inside writes all the shapes that are completely within the boundary.");
	puts( "      Boundary clips are only the min and max of a theme boundary." );
	puts( "  - Erase boundary extent or by theme boundary." );
	puts( "      Erase is the direct opposite of the Clip function." );
	puts( "  - Change coordinate value units between meters and feet.");
	puts( "      There is no way to determine the input unit of a shape file.");
	puts( "      Skip this function if the shape file is already in the correct unit.");
	puts( "      Clip and Erase will be done before the unit is changed.");
	puts( "      A shift will be done after the unit is changed."); 
	puts( "  - Shift X and Y coordinates.\n" );
	puts( "Finally, There can only be one select or unselect in the command line.");
	puts( "         There can only be one clip or erase in the command line.");
	puts( "         There can only be one unit and only one shift in the command line.\n");
	puts( "Ex: shputils in.shp out.shp   SELECT countycode 3,5,9,13,17,27");
	puts( "    shputils in.shp out.shp   CLIP   10 10 90 90 Touch   FACTOR Meter Feet");
	puts( "    shputils in.shp out.shp   FACTOR Meter 3.0");
	puts( "    shputils in.shp out.shp   CLIP   clip.shp Boundary Touch   SHIFT 40 40");
	puts( "    shputils in.shp out.shp   SELECT co 112   CLIP clip.shp Boundary Touch\n");
	puts( "USAGE: shputils  <DescribeShape>   {ALL}");
	puts( "USAGE: shputils  <InputShape>  <OutShape|AppendShape>" );
	puts( "   { <FACTOR>       <FEET|MILES|METERS|KM> <FEET|MILES|METERS|KM|factor> }" );
	puts( "   { <SHIFT>        <xshift> <yshift> }" );
	puts( "   { <SELECT|UNSEL> <Item> <valuelist> }" );
	puts( "   { <CLIP|ERASE>   <xmin> <ymin> <xmax> <ymax> <TOUCH|INSIDE|CUT> }" );
	puts( "   { <CLIP|ERASE>   <theme>      <BOUNDARY>     <TOUCH|INSIDE|CUT> }" );
	puts( "     Note: CUT is not complete and does not create intersections.");
	puts( "           For more information read programmer comment.");
	
	/****   Clip functions for Polygon and Cut is not supported
	There are several web pages that describe methods of doing this function.
	It seem easy to impliment until you start writting code.  I don't have the
	time to add these functions but a did leave a simple cut routine in the 
	program that can be called by using CUT instead of TOUCH in the 
	CLIP or ERASE functions.  It does not add the intersection of the line and
	the clip box, so polygons could look incomplete and lines will come up short.
	
	Information about clipping lines with a box:
           http://www.csclub.uwaterloo.ca/u/mpslager/articles/sutherland/wr.html
        Information about finding the intersection of two lines:
	   http://www.whisqu.se/per/docs/math28.htm
	   
THE CODE LOOKS LIKE THIS:
 ********************************************************	  
void Intersect_Lines(float x0,float y0,float x1,float y1,
                     float x2,float y2,float x3,float y3,
                     float *xi,float *yi)
                     {
//  this function computes the intersection of the sent lines
//  and returns the intersection point, note that the function assumes
//  the lines intersect. the function can handle vertical as well
//  as horizontal lines. note the function isn't very clever, it simply
//  applies the math, but we don't need speed since this is a
//  pre-processing step
//  The Intersect_lines program came from (http://www.whisqu.se/per/docs/math28.htm)

float a1,b1,c1, // constants of linear equations 
      a2,b2,c2,
      det_inv,  // the inverse of the determinant of the coefficientmatrix
      m1,m2;    // the slopes of each line
      
// compute slopes, note the cludge for infinity, however, this will
// be close enough
if ((x1-x0)!=0)
   m1 = (y1-y0)/(x1-x0);
else
   m1 = (float)1e+10;  // close enough to infinity
   
   
if ((x3-x2)!=0) 
   m2 = (y3-y2)/(x3-x2);
else
   m2 = (float)1e+10;  // close enough to infinity
   
// compute constants
a1 = m1;
a2 = m2;
b1 = -1;
b2 = -1;
c1 = (y0-m1*x0);
c2 = (y2-m2*x2);
// compute the inverse of the determinate
det_inv = 1/(a1*b2 - a2*b1);
// use Kramers rule to compute xi and yi
*xi=((b1*c2 - b2*c1)*det_inv);
*yi=((a2*c1 - a1*c2)*det_inv);
} // end Intersect_Lines
 **********************************************************/

	exit( 1 );
    }