You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
605 lines
15 KiB
605 lines
15 KiB
/****************************************************************************** |
|
* Copyright (c) 2004, Eric G. Miller |
|
* |
|
* This code is based in part on the earlier work of Frank Warmerdam |
|
* |
|
* Permission is hereby granted, free of charge, to any person obtaining a |
|
* copy of this software and associated documentation files (the "Software"), |
|
* to deal in the Software without restriction, including without limitation |
|
* the rights to use, copy, modify, merge, publish, distribute, sublicense, |
|
* and/or sell copies of the Software, and to permit persons to whom the |
|
* Software is furnished to do so, subject to the following conditions: |
|
* |
|
* The above copyright notice and this permission notice shall be included |
|
* in all copies or substantial portions of the Software. |
|
* |
|
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS |
|
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
|
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL |
|
* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
|
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING |
|
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER |
|
* DEALINGS IN THE SOFTWARE. |
|
****************************************************************************** |
|
* shpsort |
|
* |
|
* Rewrite a shapefile sorted by a field or by the geometry. For polygons, |
|
* sort by area, for lines sort by length and do nothing for all others. |
|
* |
|
* $Log: shpsort.c,v $ |
|
* Revision 1.3 2004-07-06 21:23:17 fwarmerdam |
|
* minor const warning fix |
|
* |
|
* Revision 1.2 2004/07/06 21:20:49 fwarmerdam |
|
* major upgrade .. sort on multiple fields |
|
* |
|
* Revision 1.4 2004/06/30 18:19:53 emiller |
|
* handle POINTZ, POINTM |
|
* |
|
* Revision 1.3 2004/06/30 17:40:32 emiller |
|
* major rewrite allows sorting on multiple fields. |
|
* |
|
* Revision 1.2 2004/06/23 23:19:58 emiller |
|
* use tuple copy, misc changes |
|
* |
|
* Revision 1.1 2004/06/23 21:38:17 emiller |
|
* Initial revision |
|
* |
|
* |
|
*/ |
|
|
|
#include <stdio.h> |
|
#include <stdlib.h> |
|
#include <ctype.h> |
|
#include <string.h> |
|
#include <math.h> |
|
#include "shapefil.h" |
|
|
|
enum FieldOrderEnum {DESCENDING, ASCENDING}; |
|
enum FieldTypeEnum { |
|
FIDType = -2, |
|
SHPType = -1, |
|
StringType = FTString, |
|
LogicalType = FTLogical, |
|
IntegerType = FTInteger, |
|
DoubleType = FTDouble |
|
}; |
|
|
|
struct DataUnion { |
|
int null; |
|
union { |
|
int i; |
|
double d; |
|
char *s; |
|
} u; |
|
}; |
|
|
|
struct DataStruct { |
|
int record; |
|
struct DataUnion *value; |
|
}; |
|
|
|
/* |
|
globals used in sorting, each element could have a pointer to |
|
a single data struct, but that's still nShapes pointers more |
|
memory. Alternatively, write a custom sort rather than using |
|
library qsort. |
|
*/ |
|
int nFields; |
|
int *fldIdx; |
|
int *fldOrder; |
|
int *fldType; |
|
int shpType; |
|
int nShapes; |
|
|
|
static struct DataStruct * build_index (SHPHandle shp, DBFHandle dbf); |
|
static char * dupstr (const char *); |
|
static void copy_related (const char *inName, const char *outName, |
|
const char *old_ext, const char *new_ext); |
|
static char ** split(const char *arg, const char *delim); |
|
static int compare(const void *, const void *); |
|
static double area2d_polygon (int n, double *x, double *y); |
|
static double shp_area (SHPObject *feat); |
|
static double length2d_polyline (int n, double *x, double *y); |
|
static double shp_length (SHPObject *feat); |
|
|
|
int main (int argc, char *argv[]) { |
|
|
|
SHPHandle inSHP, outSHP; |
|
DBFHandle inDBF, outDBF; |
|
int len; |
|
int i; |
|
char **fieldNames; |
|
char **strOrder = 0; |
|
struct DataStruct *index; |
|
int width; |
|
int decimals; |
|
SHPObject *feat; |
|
void *tuple; |
|
|
|
if (argc < 4) { |
|
printf("USAGE: shpsort <infile> <outfile> <field[;...]> [<(ASCENDING|DESCENDING)[;...]>]\n"); |
|
exit(EXIT_FAILURE); |
|
} |
|
|
|
inSHP = SHPOpen (argv[1], "rb"); |
|
if (!inSHP) { |
|
fputs("Couldn't open shapefile for reading!\n", stderr); |
|
exit(EXIT_FAILURE); |
|
} |
|
SHPGetInfo(inSHP, &nShapes, &shpType, NULL, NULL); |
|
|
|
/* If we can open the inSHP, open its DBF */ |
|
inDBF = DBFOpen (argv[1], "rb"); |
|
if (!inDBF) { |
|
fputs("Couldn't open dbf file for reading!\n", stderr); |
|
exit(EXIT_FAILURE); |
|
} |
|
|
|
/* Parse fields and validate existence */ |
|
fieldNames = split(argv[3], ";"); |
|
if (!fieldNames) { |
|
fputs("ERROR: parsing field names!\n", stderr); |
|
exit(EXIT_FAILURE); |
|
} |
|
for (nFields = 0; fieldNames[nFields] ; nFields++) { |
|
continue; |
|
} |
|
|
|
fldIdx = malloc(sizeof *fldIdx * nFields); |
|
if (!fldIdx) { |
|
fputs("malloc failed!\n", stderr); |
|
exit(EXIT_FAILURE); |
|
} |
|
for (i = 0; i < nFields; i++) { |
|
len = (int)strlen(fieldNames[i]); |
|
while(len > 0) { |
|
--len; |
|
fieldNames[i][len] = (char)toupper((unsigned char)fieldNames[i][len]); |
|
} |
|
fldIdx[i] = DBFGetFieldIndex(inDBF, fieldNames[i]); |
|
if (fldIdx[i] < 0) { |
|
/* try "SHAPE" */ |
|
if (strcmp(fieldNames[i], "SHAPE") == 0) { |
|
fldIdx[i] = -1; |
|
} |
|
else if (strcmp(fieldNames[i], "FID") == 0) { |
|
fldIdx[i] = -2; |
|
} |
|
else { |
|
fprintf(stderr, "ERROR: field '%s' not found!\n", fieldNames[i]); |
|
exit(EXIT_FAILURE); |
|
} |
|
} |
|
} |
|
|
|
|
|
/* set up field type array */ |
|
fldType = malloc(sizeof *fldType * nFields); |
|
if (!fldType) { |
|
fputs("malloc failed!\n", stderr); |
|
exit(EXIT_FAILURE); |
|
} |
|
for (i = 0; i < nFields; i++) { |
|
if (fldIdx[i] < 0) { |
|
fldType[i] = fldIdx[i]; |
|
} |
|
else { |
|
fldType[i] = DBFGetFieldInfo(inDBF, fldIdx[i], NULL, &width, &decimals); |
|
if (fldType[i] == FTInvalid) { |
|
fputs("Unrecognized field type in dBASE file!\n", stderr); |
|
exit(EXIT_FAILURE); |
|
} |
|
} |
|
} |
|
|
|
|
|
/* set up field order array */ |
|
fldOrder = malloc(sizeof *fldOrder * nFields); |
|
if (!fldOrder) { |
|
fputs("malloc failed!\n", stderr); |
|
exit(EXIT_FAILURE); |
|
} |
|
for (i = 0; i < nFields; i++) { |
|
/* default to ascending order */ |
|
fldOrder[i] = ASCENDING; |
|
} |
|
if (argc > 4) { |
|
strOrder = split(argv[4], ";"); |
|
if (!strOrder) { |
|
fputs("ERROR: parsing fields ordering!\n", stderr); |
|
exit(EXIT_FAILURE); |
|
} |
|
for (i = 0; i < nFields && strOrder[i]; i++) { |
|
if (strcmp(strOrder[i], "DESCENDING") == 0) { |
|
fldOrder[i] = DESCENDING; |
|
} |
|
} |
|
} |
|
|
|
/* build the index */ |
|
index = build_index (inSHP, inDBF); |
|
|
|
/* Create output shapefile */ |
|
outSHP = SHPCreate(argv[2], shpType); |
|
if (!outSHP) { |
|
fprintf(stderr, "%s:%d: couldn't create output shapefile!\n", |
|
__FILE__, __LINE__); |
|
exit(EXIT_FAILURE); |
|
} |
|
|
|
/* Create output dbf */ |
|
outDBF = DBFCloneEmpty(inDBF, argv[2]); |
|
if (!outDBF) { |
|
fprintf(stderr, "%s:%d: couldn't create output dBASE file!\n", |
|
__FILE__, __LINE__); |
|
exit(EXIT_FAILURE); |
|
} |
|
|
|
/* Copy projection file, if any */ |
|
copy_related(argv[1], argv[2], ".shp", ".prj"); |
|
|
|
/* Copy metadata file, if any */ |
|
copy_related(argv[1], argv[2], ".shp", ".shp.xml"); |
|
|
|
/* Write out sorted results */ |
|
for (i = 0; i < nShapes; i++) { |
|
feat = SHPReadObject(inSHP, index[i].record); |
|
if (SHPWriteObject(outSHP, -1, feat) < 0) { |
|
fprintf(stderr, "%s:%d: error writing shapefile!\n", __FILE__, __LINE__); |
|
exit(EXIT_FAILURE); |
|
} |
|
tuple = (void *) DBFReadTuple(inDBF, index[i].record); |
|
if (DBFWriteTuple(outDBF, i, tuple) < 0) { |
|
fprintf(stderr, "%s:%d: error writing dBASE file!\n", __FILE__, __LINE__); |
|
exit(EXIT_FAILURE); |
|
} |
|
} |
|
SHPClose(inSHP); |
|
SHPClose(outSHP); |
|
DBFClose(inDBF); |
|
DBFClose(outDBF); |
|
|
|
return EXIT_SUCCESS; |
|
|
|
} |
|
|
|
static char ** split(const char *arg, const char *delim) |
|
{ |
|
char *copy = dupstr(arg); |
|
char *cptr = copy; |
|
char **result = NULL; |
|
char **tmp; |
|
int i = 0; |
|
|
|
for (cptr = strtok(copy, delim); cptr; cptr = strtok(NULL, delim)) { |
|
tmp = realloc (result, sizeof *result * (i + 1)); |
|
if (!tmp && result) { |
|
while (i > 0) { |
|
free(result[--i]); |
|
} |
|
free(result); |
|
free(copy); |
|
return NULL; |
|
} |
|
result = tmp; |
|
result[i++] = dupstr(cptr); |
|
} |
|
|
|
free(copy); |
|
|
|
if (i) { |
|
tmp = realloc(result, sizeof *result * (i + 1)); |
|
if (!tmp) { |
|
while (i > 0) { |
|
free(result[--i]); |
|
} |
|
free(result); |
|
free(copy); |
|
return NULL; |
|
} |
|
result = tmp; |
|
result[i++] = NULL; |
|
} |
|
|
|
return result; |
|
} |
|
|
|
|
|
static void copy_related (const char *inName, const char *outName, |
|
const char *old_ext, const char *new_ext) |
|
{ |
|
char *in; |
|
char *out; |
|
FILE *inFile; |
|
FILE *outFile; |
|
int c; |
|
size_t name_len = strlen(inName); |
|
size_t old_len = strlen(old_ext); |
|
size_t new_len = strlen(new_ext); |
|
|
|
in = malloc(name_len - old_len + new_len + 1); |
|
strncpy(in, inName, (name_len - old_len)); |
|
strcpy(&in[(name_len - old_len)], new_ext); |
|
inFile = fopen(in, "rb"); |
|
if (!inFile) { |
|
free(in); |
|
return; |
|
} |
|
name_len = strlen(outName); |
|
out = malloc(name_len - old_len + new_len + 1); |
|
strncpy(out, outName, (name_len - old_len)); |
|
strcpy(&out[(name_len - old_len)], new_ext); |
|
outFile = fopen(out, "wb"); |
|
if (!out) { |
|
fprintf(stderr, "%s:%d: couldn't copy related file!\n", |
|
__FILE__, __LINE__); |
|
free(in); |
|
free(out); |
|
return; |
|
} |
|
while ((c = fgetc(inFile)) != EOF) { |
|
fputc(c, outFile); |
|
} |
|
fclose(inFile); |
|
fclose(outFile); |
|
free(in); |
|
free(out); |
|
} |
|
|
|
static char * dupstr (const char *src) |
|
{ |
|
char *dst = malloc(strlen(src) + 1); |
|
char *cptr; |
|
if (!dst) { |
|
fprintf(stderr, "%s:%d: malloc failed!\n", __FILE__, __LINE__); |
|
exit(EXIT_FAILURE); |
|
} |
|
cptr = dst; |
|
while ((*cptr++ = *src++)) |
|
; |
|
return dst; |
|
} |
|
|
|
#ifdef DEBUG |
|
static void PrintDataStruct (struct DataStruct *data) { |
|
int i, j; |
|
for (i = 0; i < nShapes; i++) { |
|
printf("data[%d] {\n", i); |
|
printf("\t.record = %d\n", data[i].record); |
|
for (j = 0; j < nFields; j++) { |
|
printf("\t.value[%d].null = %d\n", j, data[i].value[j].null); |
|
if (!data[i].value[j].null) { |
|
switch(fldType[j]) { |
|
case FIDType: |
|
case IntegerType: |
|
case LogicalType: |
|
printf("\t.value[%d].u.i = %d\n", j, data[i].value[j].u.i); |
|
break; |
|
case DoubleType: |
|
case SHPType: |
|
printf("\t.value[%d].u.d = %f\n", j, data[i].value[j].u.d); |
|
break; |
|
case StringType: |
|
printf("\t.value[%d].u.s = %s\n", j, data[i].value[j].u.s); |
|
break; |
|
} |
|
} |
|
} |
|
puts("}"); |
|
} |
|
} |
|
#endif |
|
|
|
static struct DataStruct * build_index (SHPHandle shp, DBFHandle dbf) { |
|
struct DataStruct *data; |
|
SHPObject *feat; |
|
int i; |
|
int j; |
|
|
|
/* make array */ |
|
data = malloc (sizeof *data * nShapes); |
|
if (!data) { |
|
fputs("malloc failed!\n", stderr); |
|
exit(EXIT_FAILURE); |
|
} |
|
|
|
/* populate array */ |
|
for (i = 0; i < nShapes; i++) { |
|
data[i].value = malloc(sizeof data[0].value[0] * nFields); |
|
if (0 == data[i].value) { |
|
fputs("malloc failed!\n", stderr); |
|
exit(EXIT_FAILURE); |
|
} |
|
data[i].record = i; |
|
for (j = 0; j < nFields; j++) { |
|
data[i].value[j].null = 0; |
|
switch (fldType[j]) { |
|
case FIDType: |
|
data[i].value[j].u.i = i; |
|
break; |
|
case SHPType: |
|
feat = SHPReadObject(shp, i); |
|
switch (feat->nSHPType) { |
|
case SHPT_NULL: |
|
fprintf(stderr, "Shape %d is a null feature!\n", i); |
|
data[i].value[j].null = 1; |
|
break; |
|
case SHPT_POINT: |
|
case SHPT_POINTZ: |
|
case SHPT_POINTM: |
|
case SHPT_MULTIPOINT: |
|
case SHPT_MULTIPOINTZ: |
|
case SHPT_MULTIPOINTM: |
|
case SHPT_MULTIPATCH: |
|
/* Y-sort bounds */ |
|
data[i].value[j].u.d = feat->dfYMax; |
|
break; |
|
case SHPT_ARC: |
|
case SHPT_ARCZ: |
|
case SHPT_ARCM: |
|
data[i].value[j].u.d = shp_length(feat); |
|
break; |
|
case SHPT_POLYGON: |
|
case SHPT_POLYGONZ: |
|
case SHPT_POLYGONM: |
|
data[i].value[j].u.d = shp_area(feat); |
|
break; |
|
default: |
|
fputs("Can't sort on Shapefile feature type!\n", stderr); |
|
exit(EXIT_FAILURE); |
|
} |
|
SHPDestroyObject(feat); |
|
break; |
|
case FTString: |
|
data[i].value[j].null = DBFIsAttributeNULL(dbf, i, fldIdx[j]); |
|
if (!data[i].value[j].null) { |
|
data[i].value[j].u.s = dupstr(DBFReadStringAttribute(dbf, i, fldIdx[j])); |
|
} |
|
break; |
|
case FTInteger: |
|
case FTLogical: |
|
data[i].value[j].null = DBFIsAttributeNULL(dbf, i, fldIdx[j]); |
|
if (!data[i].value[j].null) { |
|
data[i].value[j].u.i = DBFReadIntegerAttribute(dbf, i, fldIdx[j]); |
|
} |
|
break; |
|
case FTDouble: |
|
data[i].value[j].null = DBFIsAttributeNULL(dbf, i, fldIdx[j]); |
|
if (!data[i].value[j].null) { |
|
data[i].value[j].u.d = DBFReadDoubleAttribute(dbf, i, fldIdx[j]); |
|
} |
|
break; |
|
} |
|
} |
|
} |
|
|
|
#ifdef DEBUG |
|
PrintDataStruct(data); |
|
fputs("build_index: sorting array\n", stdout); |
|
#endif |
|
|
|
qsort (data, nShapes, sizeof data[0], compare); |
|
|
|
#ifdef DEBUG |
|
PrintDataStruct(data); |
|
fputs("build_index: returning array\n", stdout); |
|
#endif |
|
|
|
return data; |
|
} |
|
|
|
static int compare(const void *A, const void *B) { |
|
const struct DataStruct *a = A; |
|
const struct DataStruct *b = B; |
|
int i; |
|
int result = 0; |
|
|
|
for (i = 0; i < nFields; i++) { |
|
if (a->value[i].null && b->value[i].null) { |
|
continue; |
|
} |
|
if (a->value[i].null && !b->value[i].null) { |
|
return (fldOrder[i]) ? 1 : -1; |
|
} |
|
if (!a->value[i].null && b->value[i].null) { |
|
return (fldOrder[i]) ? -1 : 1; |
|
} |
|
switch (fldType[i]) { |
|
case FIDType: |
|
case IntegerType: |
|
case LogicalType: |
|
if (a->value[i].u.i < b->value[i].u.i) { |
|
return (fldOrder[i]) ? -1 : 1; |
|
} |
|
if (a->value[i].u.i > b->value[i].u.i) { |
|
return (fldOrder[i]) ? 1 : -1; |
|
} |
|
break; |
|
case DoubleType: |
|
case SHPType: |
|
if (a->value[i].u.d < b->value[i].u.d) { |
|
return (fldOrder[i]) ? -1 : 1; |
|
} |
|
if (a->value[i].u.d > b->value[i].u.d) { |
|
return (fldOrder[i]) ? 1 : -1; |
|
} |
|
break; |
|
case StringType: |
|
result = strcmp(a->value[i].u.s, b->value[i].u.s); |
|
if (result) { |
|
return (fldOrder[i]) ? result : -result; |
|
} |
|
break; |
|
default: |
|
fprintf(stderr, "compare: Program Error! Unhandled field type! fldType[%d] = %d\n", i, fldType[i]); |
|
break; |
|
} |
|
} |
|
return 0; |
|
} |
|
|
|
static double area2d_polygon (int n, double *x, double *y) { |
|
double area = 0; |
|
int i; |
|
for (i = 1; i < n; i++) { |
|
area += (x[i-1] + x[i]) * (y[i] - y[i-1]); |
|
} |
|
return area / 2.0; |
|
} |
|
|
|
static double shp_area (SHPObject *feat) { |
|
double area = 0.0; |
|
if (feat->nParts == 0) { |
|
area = area2d_polygon (feat->nVertices, feat->padfX, feat->padfY); |
|
} |
|
else { |
|
int part, n; |
|
for (part = 0; part < feat->nParts; part++) { |
|
if (part < feat->nParts - 1) { |
|
n = feat->panPartStart[part+1] - feat->panPartStart[part]; |
|
} |
|
else { |
|
n = feat->nVertices - feat->panPartStart[part]; |
|
} |
|
area += area2d_polygon (n, &(feat->padfX[feat->panPartStart[part]]), |
|
&(feat->padfY[feat->panPartStart[part]])); |
|
} |
|
} |
|
/* our area function computes in opposite direction */ |
|
return -area; |
|
} |
|
|
|
static double length2d_polyline (int n, double *x, double *y) { |
|
double length = 0.0; |
|
int i; |
|
for (i = 1; i < n; i++) { |
|
length += sqrt((x[i] - x[i-1])*(x[i] - x[i-1]) |
|
+ |
|
(y[i] - y[i-1])*(y[i] - y[i-1])); |
|
} |
|
return length; |
|
} |
|
|
|
static double shp_length (SHPObject *feat) { |
|
double length = 0.0; |
|
if (feat->nParts == 0) { |
|
length = length2d_polyline(feat->nVertices, feat->padfX, feat->padfY); |
|
} |
|
else { |
|
int part, n; |
|
for (part = 0; part < feat->nParts; part++) { |
|
if (part < feat->nParts - 1) { |
|
n = feat->panPartStart[part+1] - feat->panPartStart[part]; |
|
} |
|
else { |
|
n = feat->nVertices - feat->panPartStart[part]; |
|
} |
|
length += length2d_polyline (n, |
|
&(feat->padfX[feat->panPartStart[part]]), |
|
&(feat->padfY[feat->panPartStart[part]])); |
|
} |
|
} |
|
return length; |
|
} |
|
|
|
|