## Copyright (C) 2014-2018 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 ## the Free Software Foundation; either version 3 of the License, or ## (at your option) any later version. ## ## This program is distributed in the hope that it will be useful, ## but WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with this program. If not, see . ## -*- texinfo -*- ## @deftypefn {Function File} {@var{status} =} shapewrite (@var{shpstr}, @var{fname}) ## Write contents of map- or geostruct to a GIS shape file. ## ## @var{shpstr} must be a valid mapstruct or geostruct, a struct array with an ## entry for each shape feature, with fields Geometry, BoundingBox, and X and Y ## (mapstruct) or Lat and Lon (geostruct). For geostructs, Lat and Lon field ## data will be written as X and Y data. Field Geometry can have data values ## of "Point", "MultiPoint", "Line" or "Polygon", all case-insensitive. ## For each shape feature, field BoundingBox should contain the minimum and ## maximum (X,Y) coordinates in a 2x2 array [minX, minY; maxX, maxY]. ## The X and Y fields should contain X (or Latitude) and Y (or Longitude) ## coordinates for each point or vertex as row vectors; for ## poly(lines) and polygons vertices of each subfeature (if present) should be ## separated by NaN entries. ## ## M or M types (e.g., PointM, PolygonZ) can also be ## written; shapewrite.m just checks if fields "M" and/or "Z" are present in ## input mapstruct. ## ## @var{fname} should be a valid shape file name, optionally with a '.shp' ## suffix. ## ## Optional third input argument @var{atts} is one attribute name or a cellstr ## array containing a list of attribute names. Only those attribute names will ## be written to the .dbf file. Alternatively a struct can be supplied with ## attibute names contained in field "fieldName". If this argument is omitted ## all attributes will be written to the shapefile. ## ## shapewrite produces 3 files, i.e. a .shp file (the actual shape file), ## a .shx file (index file), and a .dbf file (dBase type 3) with the contents ## of additional attribute fields other than Geometry, X/Y or Lat/Lon, M, Z, ## and BoundingBox. If no additional are present, a .dbf file is created with ## the shape feature numbers as contents of field "ID". ## ## @var{status} is 1 if the complete shape file set was written successfully, ## 0 otherwise. ## ## Ref: http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf ## ## @seealso{shapedraw, shapeinfo, shaperead} ## @end deftypefn ## Author: Philip Nienhuis ## Created: 2014-12-30 function [status] = shapewrite (shp, fname="", atts=[]) status = 0; ## Input validation if (nargin < 1) print_usage; endif ## Assess shape variable type (oct or ml/geo ml/map) if (! isstruct (shp)) error ("shapewrite: [map-, geo-] struct expected for argument #2"); else ## Yep. Find out what type fldn = fieldnames (shp); if (ismember ("vals", fldn) && ismember ("shpbox", fldn)) ## Assume it is an Octave-style struct read by shaperead otype = 0; warning ("shapewrite: only Matlab-type map/geostructs can be written\n"); return; elseif (ismember ("Geometry", fldn) && all (ismember ({"X", "Y"}, fldn))) ## Assume it is a Matlab-style mapstruct otype = 1; elseif (ismember ("Geometry", fldn) && all (ismember ({"Lat", "Lon"}, fldn))) ## Assume it is a Matlab-style geostruct otype = 2; else ## Not a supported struct type error ("shapewrite: unsupported struct type.\n") endif endif ## FIXME struct field type validation if (strcmpi (shp(1).Geometry, "MultiPatch")) error ("shapewrite.m: MultiPatch type not supported"); endif ## Check for dbfwrite function if (isempty (which ("dbfwrite"))) error ("shapewrite.m: dbfwrite function not found. (io package installed \ and loaded?)"); return; endif ## Check file name if (isempty (fname)) error ("shapewrite: filename expected for input argument #2"); else [pth, fnm, ext] = fileparts (fname); if (isempty (ext)) bname = fname; fname = [fname ".shp"]; ## Later on bname.shx and bname.dbf will be created elseif (isempty (pth)) bname = fnm; else bname = [pth filesep fnm]; endif endif ## Check optional 3rd argument if (nargin > 2) if (isstruct (atts)) warning ("shapewrite.m: DbfSpec not implemented; including attributes"); ## Get attribute names from field "FieldName"; allow lowercase and camelcase f_atts = find (ismember ("fieldname", lower (fieldnames (shp)))); if (! isempty (f_atts)) atts = {shp.(fieldnames (shp (f_atts)))}; endif elseif (! iscellstr (atts) && ! ischar (atts)) error ("shapewrite.m: arg.#3: attribute name or cellstr array of attribute names expected"); endif endif ## Prepare a few things numfeat = numel (shp); if (abs (otype) >= 1) stype = strmatch (lower (shp(1).Geometry), ... {"point", "multipoint", "line", "polygon"}); if (isempty (stype)) ## Not a supported struct type error ("shapewrite: unsupported struct type.\n") else stype = [1, 8, 3, 5](stype); endif ## Preprocess geostructs if (abs (otype) == 2) ## Change Lat/Lon fields into X/Y [shp.X] = deal (shp.Lon); [shp.Y] = deal (shp.Lat); endif ## "Point" need not have a BoundingBox field => add a dummy if not found if (stype == 1 && ! ismember ("BoundingBox", fldn)) [shp.BoundingBox] = deal (repmat (zeros (2, 2), numel (shp))); endif if (any (ismember ({"M", "Z"}, fldn))) if (ismember ({"Z"}, fldn)) stype += 10; else stype += 20; endif otype = -otype; endif endif ## Only now (after input checks) open .shp and .shx files & rewind just to be sure fids = fopen (fname, "w"); if (fids < 0) error ("shapewrite: shapefile %s can't be opened for writing\n", fname); endif fseek (fids, 0, "bof"); fidx = fopen ([ bname ".shx" ], "w"); if (fidx < 0) error ("shapewrite: index file %s can't be opened for writing\n", fname); endif fseek (fidx, 0, "bof"); ## Write headers in .shp & .shx (identical). First magic number 9994 + six ## zeros, the last zero in .shp is a placeholder for the yet unknown .shp ## file length. fwrite (fids, [9994 0 0 0 0 0 0], "int32", 0, "ieee-be"); fwrite (fidx, [9994 0 0 0 0 0], "int32", 0, "ieee-be"); ## For .shx the file length in 16-bit words (single) is known: fwrite (fidx, ((numfeat * 4) + 50), "int32", 0, "ieee-be"); ## Next, shp file version fwrite (fids, 1000, "int32"); fwrite (fidx, 1000, "int32"); ## Shape feature type fwrite (fids, stype, "int32"); fwrite (fidx, stype, "int32"); ## Bounding box. Can be run later for ML type shape structs. Fill with zeros fwrite (fids, [0 0 0 0 0 0 0 0], "double"); fwrite (fidx, [0 0 0 0 0 0 0 0], "double"); ## Prepare BoundingBox limits xMin = yMin = zMin = mMin = Inf; xMax = yMax = zMax = mMax = -Inf; ## Skip to start of first record position % fseek (fids, 100, "bof"); % fseek (fidx, 100, "bof"); ## Write shape features one by one if (abs (otype) >= 1) for ishp=1:numfeat ## Write record start pos to .shx file fwrite (fidx, ftell (fids) / 2, "int32", 0, "ieee-be"); ## Prepare multipart polygons/lines. ## Find pointers to separators idx = [ 0 find(isnan (shp(ishp).X)) ]; ## Eliminate trailing NaN rows if (isnan (shp(ishp).X)) idx(end) = []; endif ## Augment idx for later on idx = unique ([ idx (numel (shp(ishp).X)+1) ]); ## Remove NaN separators idn = find (! isfinite (shp(ishp).X)); shp(ishp).X(idn) = []; shp(ishp).Y(idn) = []; if (stype >= 10) shp(ishp).Z(idn) = []; endif if (stype >= 10) shp(ishp).M(idn) = []; ## M as need not be present (we assume a NaN value then). ## Set M to a value less than -1e39 idm = find (! isfinite (shp(ishp).M)); shp(ishp).M(idm) = -1.1e39; endif ## Nr. of vertices npt = numel (shp(ishp).X); ## Pointers to parts nptr = idx(1:end-1) .- [0:numel(idx)-2]; ## Write record contents switch (stype) case 1 ## Point reclen = 10; ## Write record index number & content length (fixed) fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); fwrite (fidx, reclen, "int32", 0, "ieee-be"); ## Write shape type fwrite (fids, stype, "int32"); ## Simply write XY cordinates fwrite (fids, [shp(ishp).X shp(ishp).Y], "double"); ## Update overall BoundingBox xMin = min (xMin, shp(ishp).X); xMax = max (xMax, shp(ishp).X); yMin = min (yMin, shp(ishp).Y); yMax = max (yMax, shp(ishp).Y); case 11 ## PointZ reclen = 18; ## Write record index number & content length (fixed) fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); fwrite (fidx, reclen, "int32", 0, "ieee-be"); ## Write shape type fwrite (fids, stype, "int32"); ## Simply write XY coordinates & Z & M value fwrite (fids, [shp(ishp).X shp(ishp).Y shp(ishp).Z shp(ishp).M], "double"); ## Update overall BoundingBox xMin = min (xMin, shp(ishp).X); xMax = max (xMax, shp(ishp).X); yMin = min (yMin, shp(ishp).Y); yMax = max (yMax, shp(ishp).Y); zMin = min (zMin, shp(ishp).Z); zMax = max (zMax, shp(ishp).Z); mMin = min (mMin, shp(ishp).M); mMax = max (mMax, shp(ishp).M); case 21 ## PointM reclen = 14; ## Write record index number & content length (fixed) fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); fwrite (fidx, reclen, "int32", 0, "ieee-be"); ## Write shape type fwrite (fids, stype, "int32"); ## Simply write XY coordinates & M-value fwrite (fids, [shp(ishp).X shp(ishp).Y shp(ishp).M], "double"); ## Update overall BoundingBox xMin = min (xMin, shp(ishp).X); xMax = max (xMax, shp(ishp).X); yMin = min (yMin, shp(ishp).Y); yMax = max (yMax, shp(ishp).Y); mMin = min (mMin, shp(ishp).M); mMax = max (mMax, shp(ishp).M); case 8 ## MultiPoint reclen = (40 + 16 * numel (shp(ishp).X)) / 2; ## Write record index number & content length (fixed) fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); fwrite (fidx, reclen, "int32", 0, "ieee-be"); ## Write shape type (+4) fwrite (fids, stype, "int32"); ## Write bounding box (+32 -> 36) fwrite (fids, [shp(ishp).BoundingBox'(:)]', "double"); ## Update overall BoundingBox xMin = min (xMin, min (shp(ishp).X)); xMax = max (xMax, max (shp(ishp).X)); yMin = min (yMin, min (shp(ishp).Y)); yMax = max (yMax, max (shp(ishp).Y)); ## Write nr of points and XY data (+4 -> 40 + Nx16) fwrite (fids, numel (shp(ishp).X), "int32"); fwrite (fids, [shp(ishp).X' shp(ishp).Y']', "double"); case 18 ## MultiPointZ reclen = (72 + 32 * numel (shp(ishp).X)) / 2; ## Write record index number & content length (fixed) fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); fwrite (fidx, reclen, "int32", 0, "ieee-be"); ## Write shape type fwrite (fids, stype, "int32"); ## Write bounding box fwrite (fids, [shp(ishp).BoundingBox'(:)]', "double"); ## Update overall BoundingBox xMin = min (xMin, min (shp(ishp).X)); xMax = max (xMax, max (shp(ishp).X)); yMin = min (yMin, min (shp(ishp).Y)); yMax = max (yMax, max (shp(ishp).Y)); zMin = min (zMin, min (shp(ishp).Z)); zMax = max (zMax, max (shp(ishp).Z)); mMin = min (mMin, min (shp(ishp).M)); mMax = max (mMax, max (shp(ishp).M)); ## Write Nr of points and XY data fwrite (fids, numel (shp(ishp).X), "int32"); fwrite (fids, [shp(ishp).X' shp(ishp).Y']', "double"); ## Write Z/M range and -data in turn fwrite (fids, [min(shp(ishp).Z) max(shp(ishp).Z) shp(ishp).Z], "double"); fwrite (fids, [min(shp(ishp).M) max(shp(ishp).M) shp(ishp).M], "double"); case 28 ## MultiPointM reclen = (56 + 24 * numel (shp(ishp).X)) / 2; ## Write record index number & content length (fixed) fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); fwrite (fidx, reclen, "int32", 0, "ieee-be"); ## Write shape type fwrite (fids, stype, "int32"); ## Write bounding box fwrite (fids, [shp(ishp).BoundingBox'(:)]', "double"); ## Update overall BoundingBox xMin = min (xMin, min (shp(ishp).X)); xMax = max (xMax, max (shp(ishp).X)); yMin = min (yMin, min (shp(ishp).Y)); yMax = max (yMax, max (shp(ishp).Y)); mMin = min (mMin, min (shp(ishp).M)); mMax = max (mMax, max (shp(ishp).M)); ## Write Nr of points and XY data fwrite (fids, numel (shp(ishp).X), "int32"); fwrite (fids, [shp(ishp).X' shp(ishp).Y']', "double"); ## Write M range and M data fwrite (fids, [min(shp(ishp).M) max(shp(ishp).M) shp(ishp).M], "double"); case {3, 5} ## Polyline/-gon ## Content length reclen = (44 + (numel(idx)-1)*4 + 2*8*npt) / 2; ## Write record index number & content length (fixed) fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); fwrite (fidx, reclen, "int32", 0, "ieee-be"); ## Write shape type fwrite (fids, stype, "int32"); ## Write bounding box fwrite (fids, [shp(ishp).BoundingBox'(:)]', "double"); ## Update overall BoundingBox xMin = min (xMin, min (shp(ishp).X)); xMax = max (xMax, max (shp(ishp).X)); yMin = min (yMin, min (shp(ishp).Y)); yMax = max (yMax, max (shp(ishp).Y)); ## Write number of parts, number of points, part pointers fwrite (fids, [(numel(idx)-1) npt nptr ], "int32"); fwrite (fids, [shp(ishp).X' shp(ishp).Y']'(:), "double"); case {13, 15} ## Polyline/-gonZ ## Content length reclen = (76 + (numel(idx)-1)*4 + 32*npt) / 2; ## Write record index number & content length (fixed) fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); fwrite (fidx, reclen, "int32", 0, "ieee-be"); ## Shape type fwrite (fids, stype, "int32"); ## Write bounding box fwrite (fids, [shp(ishp).BoundingBox'(:)]', "double"); ## Update overall BoundingBox xMin = min (xMin, min (shp(ishp).X)); xMax = max (xMax, max (shp(ishp).X)); yMin = min (yMin, min (shp(ishp).Y)); yMax = max (yMax, max (shp(ishp).Y)); zMin = min (zMin, min (shp(ishp).Z)); zMax = max (zMax, max (shp(ishp).Z)); mMin = min (mMin, min (shp(ishp).M)); mMax = max (mMax, max (shp(ishp).M)); ## Write number of parts, number of points, part pointers fwrite (fids, [(numel(idx)-1) npt nptr ], "int32"); ## Write XY data fwrite (fids, [shp(ishp).X' shp(ishp).Y']'(:), "double"); fwrite (fids, [min(shp(ishp).Z) max(shp(ishp).Z) ... shp(ishp).Z], "double"); fwrite (fids, [min(shp(ishp).M) max(shp(ishp).M) ... shp(ishp).M], "double"); case {23, 25} ## Polyline/-gonM ## Content length reclen = (60 + (numel(idx)-1)*4 + 24*npt) / 2; ## Write record index number & content length (fixed) fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); fwrite (fidx, reclen, "int32", 0, "ieee-be"); ## Write shape type fwrite (fids, stype, "int32"); ## Write bounding box fwrite (fids, [shp(ishp).BoundingBox'(:)]', "double"); ## Update overall BoundingBox xMin = min (xMin, min (shp(ishp).X)); xMax = max (xMax, max (shp(ishp).X)); yMin = min (yMin, min (shp(ishp).Y)); yMax = max (yMax, max (shp(ishp).Y)); mMin = min (mMin, min (shp(ishp).M)); mMax = max (mMax, max (shp(ishp).M)); ## Write number of parts, number of points, part pointers fwrite (fids, [(numel(idx)-1) npt nptr ], "int32"); ## Write XY data fwrite (fids, [shp(ishp).X' shp(ishp).Y']'(:), "double"); ## Write M range and M data fwrite (fids, [min(shp(ishp).M) max(shp(ishp).M) ... shp(ishp).M], "double"); otherwise ## Future shape types or types unsupported yet (MultiPatch) endswitch endfor endif ## Write file length and overall BoundingBox into .shp header flen = ftell (fids); fseek (fids, 24, "bof"); fwrite (fids, flen/2, "int32", 0, "ieee-be"); fseek (fids, 36, "bof"); fwrite (fids, [xMin yMin xMax yMax], "double"); if (stype > 10) ## +-Inf & NaN not allowed in shapefiles zm = [zMin zMax mMin mMax]; zm (! isfinite (zm)) = -1e-39; fwrite (fids, zm, "double"); endif ## Close files fclose (fids); fclose (fidx); ## Write .dbf file. ## Remove basic attributes if (abs (otype) == 1) ## Attributes + shp data in mapstruct shp = rmfield (shp, {"Geometry", "BoundingBox", "X", "Y"}); elseif (abs (otype) == 2) ## Attributes + shp data in geostruct shp = rmfield (shp, {"Geometry", "BoundingBox", "Lat", "Lon", "X", "Y"}); endif if (otype < 1) shp = rmfield (shp, {"M"}); if (isfield (shp, "Z")) shp = rmfield (shp, {"Z"}); endif endif ## Write rest of attributes if (nargin == 3) ## Only write user-specified attribute selection fldn = fieldnames (shp); ## First remove regular attributes (with a value for each vertex) fldn (ismember (fldn, {"Geometry", "BoundingBox", "Lat", "Lon", "X", "Y", "M", "Z"})) = []; ## Next, only retain user-specified attributes shp = rmfield (shp, atts); endif attribs = cell (numfeat + 1, numel (fieldnames (shp))); if (! isempty (attribs)) attribs(1, :) = fieldnames (shp); attribs(2:end, :) = (squeeze (struct2cell (shp)))'; else ## Substitute ID attribute attribs{1, 1} = "ID"; [attribs{2:end}] = deal (num2cell ([1:size(shp, 2)]){:}); endif try status = dbfwrite ([ bname ".dbf"], attribs); status = 1; catch warning ("shapewrite: writing attributes to file %s failed\n", [bname ".dbf"]); end_try_catch endfunction ## Test various geometries: (1) Point %!test %! shp.Geometry = "Point"; %! shp.X = 10; %! shp.Y = 20; %! shp.Z = 30; %! shp.M = -1; %! shp.attr_1 = "Attribute1"; %! shp.attr_Z = "AttributeA"; %! shp(2).Geometry = "Point"; %! shp(2).X = 11; %! shp(2).Y = 25; %! shp(2).Z = 32; %! shp(2).M = -2; %! shp(2).attr_1 = "Attribute2"; %! shp(2).attr_Z = "AttributeB"; %! fname = tmpnam (); %! sts = shapewrite (shp, fname); %! assert (sts, 1, eps); %! ## Check index file %%! fx = fopen ([fname ".shx"], "r"); %%! fseek (fx, 100, "bof"); %%! shxinfo = fread (fx, "*int32", 0, "ieee-be"); %%! assert (shxinfo, int32 ([50; 36; 72; 36])); %%! fclose (fx); %! ## Check on filesizes, based on Esri shapewrite doc %! assert (stat ([fname ".shp"]).size, 188, 1e-10); %! assert (stat ([fname ".shx"]).size, 116, 1e-10); %! shp2 = shaperead ([fname ".shp"]); %! assert (size (shp2), [2 1]); %! flds = fieldnames (shp2); %! fields = {"Geometry", "X", "Y", "attr_1", "attr_Z"}; %! ism = ismember (fields, flds); %! ## Do we have only those fields? %! assert (numel (ism), numel (fields)); %! ## Do we have only those fields? %! assert (sum (ism), numel (ism)); %! unlink ([fname ".shp"]); %! unlink ([fname ".shx"]); %! assert ([shp2.X shp2.Y], [10 11 20 25], 1e-10); %! assert ({shp2.Geometry}, {"Point", "Point"}); %! assert ({shp2.attr_1}, {"Attribute1", "Attribute2"}); ## Test various geometries: (2) Line & Polygon %!test %! shp.Geometry = "Line"; %! shp.BoundingBox = [9 110; 19 120]; %! shp.X = [10 110 NaN 9 109]; %! shp.Y = [20 120 NaN 19 119]; %! shp.Z = [30 130 NaN 29 129]; %! shp.M = [-1 1 NaN -11 11]; %! shp.attr_1 = "Attribute1"; %! shp.attr_Z = "AttributeA"; %! shp(2).Geometry = "Line"; %! shp(2).BoundingBox = [11 211; 24 225]; %! shp(2).X = [11 111 211 NaN 11 110 NaN 12 200]; %! shp(2).Y = [25 125 225 NaN 24 124 NaN 26 200]; %! shp(2).Z = [32 132 232 NaN 31 131 NaN 33 200]; %! shp(2).M = [-2 NaN -3 NaN -22 22 NaN -33 33]; %! shp(2).attr_1 = "Attribute2"; %! shp(2).attr_Z = "AttributeB"; %! fname = tmpnam (); %! sts = shapewrite (shp, fname); %! assert (sts, 1, eps); %! ## Check index file %%! fx = fopen ([fname ".shx"], "r"); %%! fseek (fx, 100, "bof"); %%! shxinfo = fread (fx, "*int32", 0, "ieee-be"); %%! assert (shxinfo, int32 ([50; 106; 160; 156])); %%! fclose (fx); %! ## Check on filesizes, based on Esri shapewrite doc %! assert (stat ([fname ".shp"]).size, 640, 1e-10); %! assert (stat ([fname ".shx"]).size, 116, 1e-10); %! %! ## Check Matlab-style reading %! shp2 = shaperead ([fname ".shp"]); %! assert (size (shp2), [2 1]); %! flds = fieldnames (shp2); %! fields = {"Geometry", "BoundingBox", "X", "Y", "attr_1", "attr_Z"}; %! ism = ismember (fields, flds); %! ## Do we have only those fields? %! assert (numel (ism), numel (fields)); %! ## Do we have only those fields? %! assert (sum (ism), numel (ism)); %! assert ([shp2.BoundingBox], [9 110 11 211; 19 120 24 225], 1e-10); %! assert ({shp2.attr_Z}, {"AttributeA", "AttributeB"}); %! %! ## Re-read using M & Z values, read only 2nd feature %! shp2 = shaperead ([fname ".shp"], "e", "rec", 2); %! assert (size (shp2), [1 1]); %! flds = fieldnames (shp2); %! fields = {"Geometry", "BoundingBox", "X", "Y", "Z", "M", "attr_1", "attr_Z"}; %! ism = ismember (fields, flds); %! ## Do we have only those fields? %! assert (numel (ism), numel (fields)); %! ## Do we have only those fields? %! assert (sum (ism), numel (ism)); %! ## Check regular numerical values %! assert ([shp2.X; shp2.Y; shp2.Z], ... %! [11, 111, 211, NaN, 11, 110, NaN, 12, 200; ... %! 25, 125, 225, NaN, 24, 124, NaN, 26, 200; ... %! 32, 132, 232, NaN, 31, 131, NaN, 33, 200], 1e-10); %! ## Check missing M value %! assert (shp2.M, [-2, NaN, -3, NaN, -22, 22, NaN, -33, 33], 1e-10); %! assert ({shp2.attr_1, shp2.attr_Z}, {"Attribute2", "AttributeB"}); %! %! unlink ([fname ".shp"]); %! unlink ([fname ".shx"]);