summaryrefslogtreecommitdiff
path: root/inst/shaperead.m
diff options
context:
space:
mode:
Diffstat (limited to 'inst/shaperead.m')
-rw-r--r--inst/shaperead.m498
1 files changed, 275 insertions, 223 deletions
diff --git a/inst/shaperead.m b/inst/shaperead.m
index 7968730..39325b0 100644
--- a/inst/shaperead.m
+++ b/inst/shaperead.m
@@ -1,4 +1,4 @@
-## Copyright (C) 2014,2015 Philip Nienhuis
+## Copyright (C) 2014-2020 Philip Nienhuis
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
@@ -47,7 +47,7 @@
## @itemx e
## Same as 1 but M and Z type and MultiPatch shape features are accepted. The
## resulting output struct is no more ML-compatible. If the shapefile contains
-## M and/or Z type shape features the mapstruct or gestruct has extra fields M
+## M and/or Z type shape features the mapstruct or geostruct has extra fields M
## and -optionally- Z. Note that MultiPatch shape features may not have
## M-values even if Z-values are present. For MultiPatch shapes another field
## Parts is added, a Px2 array with zero-based indices to the first vertex of
@@ -97,6 +97,13 @@
## characters are significant, case doesn't matter):
##
## @table @code
+## @item Attributes
+## Normally all attributes associated with the shape features will be read
+## and returned in the output struct(s). To limit this to just some
+## attributes, enter a value consisting of a cell array of attribute names to
+## be read. To have no attributes read at all, specify @{@}, an empty cell
+## array.
+##
## @item BoundingBox
## Select only those shape items (features) whose bounding box lies within, or
## intersets in at least one point with the limits of the BoundingBox value (a
@@ -105,18 +112,20 @@
## default!
##
## @item Clip
-## (only useful in conjuction with the BoundingBox property) If a value of 1 or
-## true is supplied, clip all shapes to the bounding box limits. This option
-## may take quite a bit of processing time. If a value of "0" or false is
-## given, do not perform clipping. The default value is 0.
+## (only useful in conjuction with the BoundingBox property) If a value of 1
+## or true is supplied, clip all shapes to the bounding box limits. This
+## option may take quite a bit of processing time. If a value of "0" or false
+## is given, do not perform clipping. The default value is 0.
## Clipping is merely meant to be performed in the XY plane. Clipping 3D
-## shapes may work to some extent but can lead to unpredictable results; this
-## is especially true for MultiPatch shape types.
-## For M and Z type polylines and polygons, the M and Z values are linearly
-## interpolated for segments crossing the bounding box. As no M and Z values
-## can be computed for "new" corner nodes, NaN values are inserted there.
-## For clipping polylines and polygons the Octave-Forge geometry and octclip
-## packages need to be installed and loaded.
+## shapes is supported but may lead to unexpected results.
+## For Z and M type polylines and polygons including MultiPatch and
+## Longitude/Latitude/Height types, Z (Height) and M values for each vertex in
+## the clipped shape feature are simply copied over from the nearest vertex in
+## the original shape feature. This implies that Z and M values of new
+## vertices created on the bounding box edges may be less optimal.
+##
+## For clipping polylines and polygons the Octave-Forge geometry package needs
+## to be installed and loaded.
##
## @item Debug
## If a value of 'true' or 1 is given, shaperead echoes the current record
@@ -134,55 +143,18 @@
## (Only applicable if a Matlab-style output struct is requested). If a value
## of 'true' (or 1) is supplied, return a geostruct rather than a mapstruct.
## If a value of 0 or false is given, return a mapstruct.
-## The mere difference is that in a geostruct the fields 'X' and 'Y' are
-## replaced by 'Long' and 'Lat'. The default value is 'false' (return a
-## mapstruct').
+## The mere difference is that in a geostruct the fields 'X' and 'Y' (and
+## optionally 'Z') are replaced by 'Long' and 'Lat' (and 'Hght'). The
+## default value is 'false' (return a mapstruct').
## @end table
##
+## Ref: http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf
+##
## @seealso{geoshow, mapshow, shapedraw, shapeinfo}
## @end deftypefn
## Author: Philip Nienhuis <prnienhuis@users.sf.net>
## Created: 2014-11-07
-## Updates: **** FIXME Please leave this until the next mapping package release ***
-## 2014-11-07/10 First draft
-## 2014-11-12 Rewrite into more logical struct array
-## 2014-11-13 Fix indexing in NaN row insertion for outopts=1
-## 2014-11-14 Added Matlab-compatible output mapstruct option
-## 2014-11-15 Removed redundant code leading to a 75 X speed increase
-## '' Fix ML-style X and Y coordinates
-## 2014-11-16 Add shape type to each struct element for ML-style output
-## '' Initial input options parsing
-## '' RecordNumbers, UseGeoCoords and BoundingBox implemented
-## 2014-11-17 Double output structs (ML-compatible)
-## 2014-11-18 Polygon clipping & test for geometry package
-## 2014-11-26 Polygon clipping debugged & working
-## 2014-11-28 Return empty cell arr if no shape features read
-## '' Return full mapstruct array for no output arg
-## 2014-11-29 Doubly run inpolygon to include features > and < BoundingBox
-## 2014-11-30 Implement Polygon & Point clipping
-## 2014-12-15 Last debug of Polyline clipping
-## 2014-12-17 Fix MultiPatch reading, esp. detection of missing M (measure)
-## 2014-12-29 Implement selecting cellstr attributes in .dbf file
-## 2014-12-30 Do not transpose npr values for MultiPatch
-## 2015-01-01 Catch null shapes and prevent reading their attributes
-## 2015-01-10 Use .shx (if present) to speed up RecordNumbers searches
-## 2015-01-12 For outopts=1, apply cell2mat to numeric & logical attributes
-## 2015-01-14 Swap outopts values for ml (now 0) and plt style (now 3). In
-## passing, fix wrong switch case id at ~L.651.
-## '' Reserve outopts = 1 for future use (ML-comp. mapstruct w. Z & M
-## atts & multipatch feature type
-## 2015-01-19 Do not include/process NULL shape types once record has been read
-## 2015-01-27 Stricter checks on ML compatibility for outopts=1; fix some
-## wrong references to outopts type, add ext/e for outopts=1;
-## update texinfo header; speed-up reading MultiPatch
-## 2015-04-19 Make warning setting local
-## '' Add filesep to basename
-## 2015-07-10 Separate checks for octclip and geometry package
-## '' Make rest of warning settings local
-## '' Avoid leading filesep in case of no full path name in basename
-## '' Improve struct type check and error message
-## **** FIXME Please leave this until the next mapping package release ***
function [ outs, oatt ] = shaperead (fname, varargin);
@@ -215,59 +187,75 @@ function [ outs, oatt ] = shaperead (fname, varargin);
outopts = 0;
elseif (nargin == 2)
## Assume filename + outopts was supplied
- outopts = varargin{1};
+ if (isempty (varargin{1}))
+ ## Assume ML-style output
+ outopts = 0;
+ else
+ outopts = varargin{1};
+ endif
varargin = {};
elseif (rem (nargin, 2) == 0)
## Even number of input args => Outopts & at least one pair of varargin
outopts = varargin{1};
varargin(1) = [];
else
- ## Only filename and varargin supplied
- outopts = 0;
- ## Just to be sure check, if output type was selected anyway; it would
- ## imply a missing opts value
- if (any (strncmp (varargin{1}, {"ml", "ext", "oct", "dat"}, 1)) || ...
- isnumeric (varargin{1}))
- outopts = varargin{1};
- varargin(1) = [];
+ ## Odd nr; maybe only filename and prop/val(s) supplied, outstyle skipped
+ if (ischar (varargin{1}))
+ ## Check arg#2, must be a property name then
+ if (! ismember (lower (varargin{1}(1:min(3, numel (varargin{1})))), ...
+ {"att", "bou", "cli", "deb", "rec", "sel", "use"}))
+ error ("shaperead: property name expected for arg. #2");
+ endif
+ else
+ ## no outstyle or property => wrong input
+ print_usage ();
endif
+ outopts = 0;
endif
## Check output type arg
if (isnumeric (outopts))
if (outopts < 0 || outopts > 3)
- error ("shaperead: value for arg. #2 out of range 0-3\n");
+ error ("shaperead: arg. #2 integer value out of range 0-3\n");
endif
elseif (ischar (outopts))
outopts = lower (outopts);
if (! any (strncmp (outopts, {"ml", "ext", "oct", "dat"}, 1)))
- error ("shaperead: value for arg. #2 should be one of 'ml', 'ext', 'oct' or 'dat'\n");
+ error ("shaperead: arg. #2 char value should be one of 'ml', 'ext', \
+'oct' or 'dat'\n");
endif
switch outopts
case {"ml", "m"}
outopts = 0;
case {"ext", "e"}
- outopts = 0;
+ outopts = 1;
case {"oct", "o"}
outopts = 2;
case {"dat", "d"}
outopts = 3;
otherwise
- error ("shaperead: illegal value for arg. #2: '%s' - expected 'ml', 'ext', 'oct' of 'dat'", ...
+ error ("shaperead: illegal value for arg. #2: '%s' - expected 'ml', \
+'ext', 'oct' of 'dat'", ...
outopts);
endswitch
else
error ("shaperead: numeric or character type expected for arg. #2\n");
endif
- ## ---------------------- 1. .shx ----------------------
-
+ ## Check .shp file existence
+ fidp = fopen (fname, "r");
+ ## Postpone file opening error until after other provisional input validation
+ if (fidp > 0)
+ ## Temporarily close to avoid file handle leaks during further input checks
+ fclose (fidp);
+ endif
+ ## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Open .shx file to help speed up seeks to next records. We need this info
## for s_recs check below
have_shx = 0;
fidx = fopen ([bname ".shx"], "r");
if (fidx < 0)
- warning ("shaperead: index file %s not found\n", [fnm ".shx"]);
+ s_recs = [];
else
fseek (fidx, 24, "bof");
fxlng = fread (fidx, 1, "int32", 72, "ieee-be");
@@ -278,79 +266,86 @@ function [ outs, oatt ] = shaperead (fname, varargin);
fclose (fidx);
## Get indices & lengths in bytes
ridx *= 2;
+ have_shx = 1;
+ s_recs = [1 : size(ridx, 1)];
endif
## Parse options; first set defaults
clip = 0;
dbug = 0;
- s_atts = {};
+ s_atts = [];
s_bbox = [];
s_geo = 0;
- s_recs = [];
## Init collection of records that meet BB criteria
bb_union = [];
## Process input args
for ii = 1:2:length (varargin)
- switch (lower (varargin{ii})(1:3))
- case "att"
- ## Select records based on attribute values
- s_atts = varargin{ii+1};
- case "bou"
- ## Select whether record / shape feature points lie inside or outside limits
- try
- s_bbox = double (varargin{ii+1});
- if (numel (s_bbox) != 4)
- error ("shaperead: 2 X 2 numeric array expected for BoundingBox\n");
- endif
- catch
- error ("shaperead: numeric2 X 2 array expected for BoundingBox\n");
- end_try_catch
- ## Initialize supplementary polygon array here
- sbox = [s_bbox(1) s_bbox(2); s_bbox(1) s_bbox(4); s_bbox(3) s_bbox(4); ...
- s_bbox(3) s_bbox(2); s_bbox(1) s_bbox(2)];
- case "cli"
- ## Clip polygons to requested BoundingBox
- try
- clip = logical (varargin{ii+1});
- catch
- error ("numeric or logical value expected for 'Clip'\n");
- end_try_catch
- case "deb"
- ## Set verbose output (count records)
- try
- dbug = logical (varargin{ii+1});
- catch
- error ("numeric or logical value expected for 'Debug'\n");
- end_try_catch
- case "rec"
- ## Select record nrs directly. Check for proper type & clean up
- try
- s_recs = sort (unique (double (varargin{ii+1}(:))));
- have_shx = fidx > 0;
- if (have_shx && any (s_recs > nrec))
- printf ("shaperead: requ. record nos. > nr. of records (%d) ignored\n", ...
- nrec);
- s_recs (find (srecs > nrec)) = [];
+ if (! ischar (varargin{ii}))
+ error ("shaperead: property %d: property name expected but got a %s value", ...
+ (ii+1)/2, class (varargin{ii}));
+ elseif (numel (varargin{ii}) < 3)
+ warning ("shaperead: unknown option '%s' - ignored\n", varargin{ii});
+ else
+ switch (lower (varargin{ii})(1:3))
+ case "att"
+ ## Select records based on attribute values
+ s_atts = varargin{ii+1};
+ case "bou"
+ ## Select whether record/shape features partly lie inside or outside limits
+ try
+ s_bbox = double (varargin{ii+1});
+ if (numel (s_bbox) != 4)
+ error ("shaperead: 2 X 2 numeric array expected for BoundingBox\n");
+ endif
+ catch
+ error ("shaperead: numeric 2 X 2 array expected for BoundingBox\n");
+ end_try_catch
+ ## Initialize supplementary polygon array here
+ sbox = [s_bbox(1) s_bbox(3); s_bbox(2) s_bbox(3); s_bbox(2) s_bbox(4); ...
+ s_bbox(1) s_bbox(4); s_bbox(1) s_bbox(3)];
+ case "cli"
+ ## Clip polygons to requested BoundingBox
+ try
+ clip = logical (varargin{ii+1});
+ catch
+ error ("numeric or logical value expected for 'Clip'\n");
+ end_try_catch
+ case "deb"
+ ## Set verbose output (count records)
+ try
+ dbug = logical (varargin{ii+1});
+ catch
+ error ("numeric or logical value expected for 'Debug'\n");
+ end_try_catch
+ case "rec"
+ ## Select record nrs directly. Check for proper type & clean up
+ try
+ s_recs = sort (unique (double (varargin{ii+1}(:))));
+ if (have_shx && any (s_recs > nrec))
+ printf ("shaperead: requ. record nos. > nr. of records (%d) ignored\n", ...
+ nrec);
+ s_recs (find (srecs > nrec)) = [];
+ endif
+ catch
+ error ("shaperead: numeric value or array expected for RecordNumbers\n");
+ end_try_catch
+ case "sel"
+ ## A hard one, to be implemented later?
+ printf ("shaperead: 'Selector' option not implemented, option ignored\n");
+ case "use"
+ ## Return a geostruct or a mapstruct (default). Only for ML-structs
+ if (outopts != 0)
+ error ("shaperead: UseGeoCoords only valid for Matlab-style output\n");
endif
- catch
- error ("shaperead: numeric value or array expected for RecordNumbers\n");
- end_try_catch
- case "sel"
- ## A hard one, to be implemented later?
- printf ("shaperead: 'Selector' option not implemented, option ignored\n");
- case "use"
- ## Return a geostruct or a mapstruct (default). Only for ML-structs
- if (outopts != 0)
- error ("shaperead: UseGeoCoords only valid for Matlab-style output\n");
- endif
- try
- s_geo = logical (varargin{ii+1});
- catch
- error ("shaperead: logical value type expected for 'UseGeoCoords'\n");
- end_try_catch
- otherwise
- warning ("shaperead: unknown option '%s' - ignored\n", varargin{ii});
- endswitch
+ try
+ s_geo = logical (varargin{ii+1});
+ catch
+ error ("shaperead: logical value type expected for 'UseGeoCoords'\n");
+ end_try_catch
+ otherwise
+ warning ("shaperead: unknown option '%s' - ignored\n", varargin{ii});
+ endswitch
+ endif
endfor
## Post-processing
if (clip)
@@ -358,10 +353,10 @@ function [ outs, oatt ] = shaperead (fname, varargin);
warning ("shaperead: no BoundingBox supplied => Clip option ignored.\n");
clip = 0;
endif
- if (isempty (which ("oc_polybool")))
- ## No OF octclip package?
- printf ("shaperead: function 'oc_polybool' not found. Clip option ignored\n");
- warning (" (OF octclip package installed and loaded?)\n");
+ if (isempty (which ("clipPolygon_clipper")))
+ ## No OF geometry package?
+ printf ("shaperead: function 'clipPolygon' not found. Clip option ignored\n");
+ warning (" (OF geometry package installed and loaded?)\n");
clip = 0;
endif
if (isempty (which ("distancePoints")))
@@ -372,14 +367,21 @@ function [ outs, oatt ] = shaperead (fname, varargin);
endif
endif
+ if (fidp < 0)
+ ## Only now convey file open error message
+ error ("shaperead: can't open file %s\n", fname);
+ else
+ ## Open .shp file
+ fidp = fopen (fname, "r");
+ endif
+ if (fidx < 0)
+ warning ("shaperead: index file %s not found\n", [fnm ".shx"]);
+ endif
+
## ============= Preparations done, now we can start reading ============
## ---------------------- 2. Read .shp file proper ----------------------
## Start reading header
- fidp = fopen (fname, 'r');
- if (fidp < 0)
- error ("shaperead: couldn't open file %s\n", fname);
- endif
fseek (fidp, 0, "bof");
## Read & check file code
@@ -405,12 +407,13 @@ function [ outs, oatt ] = shaperead (fname, varargin);
## it. Initial tries showed no significant speed advantage yet over
## the incremental allocation scheme implemented in Ls. 611+ & 631+
## for oct/plt style output. For ml-style it is much harder as
- ## we'd need to preallocate a potentially very heterogeneous strtuct
+ ## we'd need to preallocate a potentially very heterogeneous struct
## array.
## Prepare for unsupported (in ML output) shape types
unsupp = 0;
- ign_mz = (outopts == 0);
+ ## Echo warning if dbug was set
+ ign_mz = (outopts == 0) && dbug;
## Buffer for record data
BUFSIZE = 10000;
## Read records, 1 by 1. Initialize final array
@@ -419,7 +422,7 @@ function [ outs, oatt ] = shaperead (fname, varargin);
nshp = 0;
## Temp pointer to keep track of array size and increase it w. BUFSIZE rows
ivals = 1;
- ## Assume file has M (measure) values
+ ## Provisionally assume file has M (measure) values
has_M = true;
## Record index number (equals struct element number)
@@ -469,40 +472,50 @@ function [ outs, oatt ] = shaperead (fname, varargin);
case 8 ## Multipoint
## Read bounding box
- tbbox = fread (fidp, 4, "double");
+ tbbox(1:4) = fread (fidp, 4, "double");
tnpt = fread (fidp, 1, "int32");
- val = reshape (fread (fidp, tnpt*2, "double"), 2, [])';
+ val(1:tnpt, 1:2) = reshape (fread (fidp, tnpt*2, "double"), 2, [])';
val(1:tnpt, 3:4) = NaN;
## Copy rec index & type down
- val(2:tnpt, 5:6) = val(1, 5:6);
+ val(2:tnpt, 5:6) = repmat (val(1, 5:6), tnpt-1, 1);
tnpr = 0;
case 18 ## MultipointZ
- tbbox = fread (fidp, 4, "double");
+ tbbox(1:4) = fread (fidp, 4, "double");
tnpt = fread (fidp, 1, "int32");
val(1:tnpt, 1:2) = reshape (fread (fidp, tnpt*2, "double"), 2, [])';
## Z min & max values
tbbox(5:6) = fread (fidp, 2, "double");
## Augment val array with Z values
- val(1:tnpt, 3) = [ rec(ir).points fread(fidp, tnpt, "double")' ];
+ val(1:tnpt, 3) = fread(fidp, tnpt, "double")';
## M min & max values
tbbox(7:8) = fread (fidp, 2, "double");
- ## Augment val array with M values
- val(1:tnpt, 4) = [ rec(ir).points fread(fidp, tnpt, "double")' ];
- ## Copy rec index & type down
- val(2:tnpt, 5:6) = repmat (val(1, 5:6), tnpt-1, 1);
- tnpr = 0;
+ if (val(1, 5) == 1)
+ has_M = checkM (fidp, val(1, 6), shpbox);
+ endif
+ if (has_M)
+ ## Augment val array with M values
+ val(1:tnpt, 4) = fread(fidp, tnpt, "double")';
+ ## Copy rec index & type down
+ val(2:tnpt, 5:6) = repmat (val(1, 5:6), tnpt-1, 1);
+ tnpr = 0;
+ endif
case 28 ## MultipointM
tbbox(1:4) = fread (fidp, 4, "double");
tnpt = fread (fidp, 1, "int32");
val(1:tnpt, 1:2) = reshape (fread (fidp, tnpt*2, "double"), 2, [])';
## Insert empty column for Z
- val(1:tnpt, 4) = NaN;
- ## M min & max values
- bbox(7:8) = fread (fidp, 2, "double");
- ## Augment val array with Z values
- val(1:tnpt, 4) = [ rec(ir).points fread(fidp, tnpt, "double")' ];
+ val(1:tnpt, 3:4) = NaN;
+ if (val(1, 5) == 1)
+ has_M = checkM (fidp, val(1, 6), shpbox);
+ endif
+ if (has_M)
+ ## M min & max values
+ tbbox(7:8) = fread (fidp, 2, "double");
+ ## Augment val array with M values
+ val(1:tnpt, 4) = fread(fidp, tnpt, "double")';
+ endif
## Copy rec index & type down
val(2:tnpt, 5:6) = repmat (val(1, 5:6), tnpt-1, 1);
tnpr = 0;
@@ -531,11 +544,16 @@ function [ outs, oatt ] = shaperead (fname, varargin);
## Z min & max values + data
tbbox(5:6) = fread (fidp, 2, "double");
val(1:tnpt, 3) = fread(fidp, tnpt, "double")';
- ## M min & max values + data
- tbbox(7:8) = fread (fidp, 2, "double");
- val(1:tnpt, 4) = fread(fidp, tnpt, "double")';
- ## Copy rec index and type down
- val(2:tnpt, 5:6) = repmat (val(1, 5:6), tnpt-1, 1);
+ if (val(1, 5) == 1)
+ has_M = checkM (fidp, val(1, 6), shpbox);
+ endif
+ if (has_M)
+ ## M min & max values + data
+ tbbox(7:8) = fread (fidp, 2, "double");
+ val(1:tnpt, 4) = fread(fidp, tnpt, "double")';
+ ## Copy rec index and type down
+ val(2:tnpt, 5:6) = repmat (val(1, 5:6), tnpt-1, 1);
+ endif
case {23, 25} ## Polyline/-gonM
tbbox(1:4) = fread (fidp, 4, "double");
@@ -547,11 +565,16 @@ function [ outs, oatt ] = shaperead (fname, varargin);
val(1:tnpt, 1:2) = reshape (fread (fidp, tnpt*2, "double"), 2, [])';
## No Z data
val(1:tnpt, 3) = NaN;
- ## M min & max values + data
- tbbox(7:8) = fread (fidp, 2, "double");
- val(1:tnpt, 4) = fread(fidp, tnpt, "double")';
- ## Copy rec index and type down
- val(2:tnpt, 5:6) = repmat (val(1, 5:6), tnpt-1, 1);
+ if (val(1, 5) == 1)
+ has_M = checkM (fidp, val(1, 6), shpbox);
+ endif
+ if (has_M)
+ ## M min & max values + data
+ tbbox(7:8) = fread (fidp, 2, "double");
+ val(1:tnpt, 4) = fread(fidp, tnpt, "double")';
+ ## Copy rec index and type down
+ val(2:tnpt, 5:6) = repmat (val(1, 5:6), tnpt-1, 1);
+ endif
case 31 ## Multipatch
tbbox(1:4) = fread (fidp, 4, "double");
@@ -562,16 +585,7 @@ function [ outs, oatt ] = shaperead (fname, varargin);
## Provisionally transpose, this is later on reset after NaN insertion
tnpr = reshape (fread (fidp, nparts*2, "int32")', [], 2)';
## Read XY point coordinates
-% val = NaN (tnpt + size (tnpr, 2) - 1, 6);
-% ipt = 1;
-% ttnpr = [tnpr(1, :) tnpt];
-% for ii=2:numel (ttnpr)
-% val(ipt:ipt+ttnpr(ii)-1, 1:2) = ...
-% reshape (fread (fidp, ttnpr(ii)*2, "double"), 2, [])';
-% ipt +=ttnpr(ii) + 1;
-% endfor
val(1:tnpt, 1:2) = reshape (fread (fidp, tnpt*2, "double"), 2, [])';
-
## Z min & max values + data. Watch out for incomplete .shp file
EOF = (ftell (fidp) > flngt-2);
if (! EOF)
@@ -581,19 +595,7 @@ function [ outs, oatt ] = shaperead (fname, varargin);
fptr = ftell (fidp);
EOF = (fptr > flngt-2);
if (! EOF && has_M)
- ## Check if the file really contains M values: try to read ("peek")
- ## three int32 record numbers matching the next record nr. to be
- ## read, the next record length and 31 (= MultiPatch type)
- tmp = fread (fidp, 2, "int32", 0, "ieee-be");
- tmp(3) = fread (fidp, 1, "int32");
- fseek (fidp, fptr, "bof");
- if (tmp(1) == (val(1, 5) + 1) && tmp(3) == 31)
- ## Undoubtedly a record index number.
- has_M = false;
- ## Empirically the min/maxM values in file header mean nothing, so
- shpbox.M(1) = 0;
- shpbox.M(2) = 0;
- endif
+ has_M = checkM (fidp, val(1, 6), shpbox);
endif
## M min & max values + data. Watch out for incomplete .shp file
if (! EOF && has_M)
@@ -615,9 +617,17 @@ function [ outs, oatt ] = shaperead (fname, varargin);
endswitch
+ ## Check if (X, Y) are valid coordinates
+ if (any (abs (val(:, 1:2)) > 1.797e308))
+ ## Probably +/- Inf
+ rincl = 0;
+ if (dbug)
+ printf ("shape# %d has no finite XY coordinates, skipped\n", ir);
+ endif
+
## Detect if shape lies (partly) within or completely out of BoundingBox.
## Null shapes are automatically skipped
- if (! isempty (s_bbox))
+ elseif (! isempty (s_bbox))
## Just check if any shape feature bounding box corner lies in s_bbox
tbox = [tbbox(1) tbbox(2); tbbox(1) tbbox(4); ...
tbbox(3) tbbox(4); tbbox(3) tbbox(2); tbbox(1) tbbox(2)];
@@ -649,14 +659,10 @@ function [ outs, oatt ] = shaperead (fname, varargin);
if (rincl && clip)
## What to do depends on shape type. Null and MultiPatch aren't clipped
switch val(1, 6)
- case {3, 13, 33} ## Polyline
- ## Temporarily silence Octave a bit, then call clipplg
+ case {3, 13, 23, 5, 15, 25, 31} ## Polyline/gon, Multipatch
+ ## Temporarily silence Octave a bit, then call clippln
warning ("off", "Octave:broadcast", "local");
[val, tnpt, tnpr] = clippln (val, tnpt, tnpr, sbox, val(1, 6));
- case {5, 15, 25, 31} ## Polygon, Multipatch
- ## Temporarily silence Octave a bit, then call clipplg
- warning ("off", "Octave:broadcast", "local");
- [val, tnpt, tnpr] = clipplg (val, tnpt, tnpr, sbox, val(1, 6));
otherwise
warning ("shaperead: unknown shape type found (%d) - ignored\n", ...
val(1, 6));
@@ -682,6 +688,10 @@ function [ outs, oatt ] = shaperead (fname, varargin);
## Keep track of nr of shape features read
++nshp;
+ ## M-values < -1e39 really mean absent values
+ im = find (val(:, 4) < -1e39);
+ val(im, 4) = NaN;
+
if (outopts < 3)
## Prepare an Octave or Matlab style struct optimized for fast
## plotting by inserting a NaN row after each polyline/-gon part
@@ -700,7 +710,7 @@ function [ outs, oatt ] = shaperead (fname, varargin);
## Shape either included by default or it lies in requested BoundingBox
switch outopts
- case {0, 1, "ml", "e"}
+ case {0, 1}
## Return a ML compatible mapstruct. Attributes will be added later
if (ign_mz && val(1, 6) >= 10)
printf ("shaperead: M and Z values ignored for ml-style output\n");
@@ -708,55 +718,60 @@ function [ outs, oatt ] = shaperead (fname, varargin);
endif
switch val(1, 6)
case {1, 11, 21} ## Point
- outs(end+1).Geometry = "Point";
+ outs(end+1, 1).Geometry = "Point";
case {8, 18, 28} ## Multipoint
- outs(end+1).Geometry = "MultiPoint";
+ outs(end+1, 1).Geometry = "MultiPoint";
case {3, 13, 23} ## Polyline
- outs(end+1).Geometry = "Line";
+ outs(end+1, 1).Geometry = "Line";
case {5, 15, 25} ## Polygon
- outs(end+1).Geometry = "Polygon";
+ outs(end+1, 1).Geometry = "Polygon";
otherwise
if (outopts == 1)
## "Extended" ML-style output struct
if (val(1, 6) == 31) ## MultiPatch
- outs(end+1).Geometry = "MultiPatch";
- outs(end).Parts = tnpr;
+ outs(end+1, 1).Geometry = "MultiPatch";
+ outs(end, 1).Parts = tnpr;
endif
else
if (! unsupp)
- warning ("shaperead: shapefile contains unsupported shape types\n");
- outs = [];
+ warning ("shaperead: shapefile contains unsupported shape \
+types\n");
+ outs = oatt = [];
return
endif
- outs(end+1).Geometry = val(1, 6);
+ outs(end+1, 1).Geometry = val(1, 6);
endif
endswitch
- outs(end).BoundingBox = reshape (tbbox(1:4), 2, [])';
+ ## Omit BoundingBox for Point
+ if (all ([1, 11, 21] - val(1, 6)))
+ outs(end).BoundingBox = reshape (tbbox(1:4), 2, [])';
+ endif
if (s_geo)
- outs(end).Lon = val(:, 1)';
- outs(end).Lat = val(:, 2)';
+ outs(end, 1).Lon = val(:, 1)';
+ outs(end, 1).Lat = val(:, 2)';
else
- outs(end).X = val(:, 1)';
- outs(end).Y = val(:, 2)';
+ outs(end, 1).X = val(:, 1)';
+ outs(end, 1).Y = val(:, 2)';
endif
## (ML-incompatible) add Z- and optional M-values, if any
if (outopts == 1 && any (isfinite (val(:, 4))))
- outs(end).M = val(:, 4)';
+ outs(end, 1).M = val(:, 4)';
## FIXME Decision needed if field Geometry should reflect the type
## outs{end}.Geometry = [ outs{end}.Geometry "M"];
endif
if (outopts == 1 && any (isfinite (val(:, 3))))
- outs(end).Z = val(:, 3)';
+ outs(end, 1).Z = val(:, 3)';
## FIXME Decision needed if field Geometry should reflect the type
## outs{end}.Geometry = [ outs{end}.Geometry "Z"];
endif
## (ML-incompatible) add Z- and optional M-values, if any
if (dbug)
## Add a field with shape feature identifier for boundingbox
- outs(end).___Shape_feature_nr___ = val(1, 5);
+ outs(end, 1).___Shape_feature_nr___ = val(1, 5);
endif
- case {2, "oct"}
+ case {2}
+ ## Return an Octave style struct.
## Append to vals array; keep track of appended nr. of rows
lvals = size (vals, 1);
lval = size (val, 1);
@@ -779,8 +794,8 @@ function [ outs, oatt ] = shaperead (fname, varargin);
endif
bbox = [bbox; tbbox];
- case {3, "dat"}
- ## Return an Octave style struct.
+ case {3}
+ ## Return a compressed Octave style struct.
## Simply append to vals array; keep track of appended nr. of rows
lvals = size (vals, 1);
lval = size (val, 1);
@@ -815,12 +830,13 @@ function [ outs, oatt ] = shaperead (fname, varargin);
until (ftell (fidp) > flngt - 3 ||
(have_shx && ir > numel (s_recs))); ## i.e., within last
- ## file bytes or
+ ## file bytes or all
+ ## req. records read
fclose (fidp);
## If no shape was read, or none fitted within BoundingBox, return empty
if (nshp < 1)
- outs = {};
+ outs = oatt = {};
return;
endif
@@ -875,11 +891,16 @@ function [ outs, oatt ] = shaperead (fname, varargin);
endif
## ---------------------- 3. .dbf ----------------------
-
+
+ if (iscell (s_atts) && isempty (s_atts))
+ ## {} indicates no attributes to be read. Morph into ""
+ return;
+ endif
## Check if dbfread is available
if (isempty (which ("dbfread")))
printf ("shaperead.m: dbfread function not found. No attributes will be added.\n");
printf (" (io package installed and loaded?)\n");
+ oatt = {};
else
## Try to read the .dbf file
@@ -887,8 +908,7 @@ function [ outs, oatt ] = shaperead (fname, varargin);
atts = dbfread ([ bname ".dbf" ], s_recs, s_atts);
if (outopts < 2)
## First check if any attributes match fieldnames; if so, pre-/append "_"
- tags = {"shpbox", "vals", "bbox", "npt", "npr", "idx", "Geometry", ...
- "BoundingBox", "X", "Y", "Lat", "Lon"};
+ tags = {"Geometry", "BoundingBox", "X", "Y", "Lat", "Lon"};
for ii=1:numel (tags)
idx = find (strcmp (tags{ii}, atts(1, :)));
if (! isempty (idx))
@@ -908,9 +928,13 @@ function [ outs, oatt ] = shaperead (fname, varargin);
for ii=1:size (atts, 2)
[oatt.(atts{1, ii})] = deal (atts(2:end, ii){:});
endfor
+ oatt = oatt';
endif
else
## Octave output struct. Add attributes columns as struct fields
+ ## Check if any attributes match fieldnames; if so, pre-/append "_"
+ tags = {"shpbox", "vals", "bbox", "npt", "npr", "idx", "Geometry", ...
+ "BoundingBox", "X", "Y", "Lat", "Lon"};
for ii=1:size (atts, 2)
outs.fields(end+1) = atts{1, ii};
if (islogical (atts{2, ii}) || isnumeric (atts{2, ii}))
@@ -928,3 +952,31 @@ function [ outs, oatt ] = shaperead (fname, varargin);
endif
endfunction
+
+
+function has_M = checkM (fidp, ft, shpbox)
+
+ has_M = true;
+ ## Check if we do have M-values. If so, the next 3 4-byte words
+ ## comprise the next record nr, rec length and shape type
+ fptr = ftell (fidp);
+ tmp = fread (fidp, 2, "int32", 0, "ieee-be");
+ tmp(3) = fread (fidp, 1, "int32");
+ if (tmp(1) == 2 && tmp(3) == ft)
+ ## Looks like next record header + next shape feature type => no M
+ has_M = false;
+ shpbox.M(1) = 0;
+ shpbox.M(2) = 0;
+ endif
+ fseek (fidp, fptr, "bof");
+
+endfunction
+
+
+## Only check input validation. I/O is tested in shapewrite.m
+%!error <arg. #2 char value should be one of> shaperead ('tst.shp', 'j');
+%!error <shaperead: arg. #2 integer value out of range> shaperead ('tst.shp', 7);
+%!error <illegal value for arg. #2:> shaperead ('tst.shp', 'deb')
+%!error < property name expected> shaperead ('tst.shp', "ml", []);
+%!error <numeric or logical value expected> shaperead ('tst.shp', 'deb', {});
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-17 01:43:19 +0000