1## Copyright (C) 2014-2020 Philip Nienhuis 2## 3## This program is free software; you can redistribute it and/or modify it 4## under the terms of the GNU General Public License as published by 5## the Free Software Foundation; either version 3 of the License, or 6## (at your option) any later version. 7## 8## This program is distributed in the hope that it will be useful, 9## but WITHOUT ANY WARRANTY; without even the implied warranty of 10## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11## GNU General Public License for more details. 12## 13## You should have received a copy of the GNU General Public License 14## along with this program. If not, see <http://www.gnu.org/licenses/>. 15 16## -*- texinfo -*- 17## @deftypefn {Function File} {@var{status} =} shapewrite (@var{shpstr}, @var{fname}) 18## @deftypefnx {Function File} {@var{status} =} shapewrite (@var{shpstr}, @var{fname}, @var{atts}) 19## Write contents of map- or geostruct to a GIS shape file. 20## 21## @var{shpstr} must be a valid mapstruct or geostruct, a struct array with an 22## entry for each shape feature, with fields Geometry, BoundingBox, and X and Y 23## (mapstruct) or Lat and Lon (geostruct). For geostructs, Lat and Lon field 24## data will be written as X and Y data. Field Geometry can have data values 25## of "Point", "MultiPoint", "Line" or "Polygon", all case-insensitive. 26## For each shape feature, field BoundingBox should contain the minimum and 27## maximum (X,Y) coordinates in a 2x2 array [minX, minY; maxX, maxY]. 28## The X and Y fields should contain X (or Latitude) and Y (or Longitude) 29## coordinates for each point or vertex as row vectors; for 30## poly(lines) and polygons vertices of each subfeature (if present) should be 31## separated by NaN entries. 32## 33## <Geometry>M or <Geometry>Z types (e.g., PointM, PolygonZ) can also be 34## written; shapewrite.m just checks if fields "M" and/or "Z" are present in 35## input mapstruct. 36## 37## @var{fname} should be a valid shape file name, optionally with a '.shp' 38## suffix. 39## 40## Optional third input argument @var{atts} is one attribute name or a cellstr 41## array containing a list of attribute names; only those attributes will be 42## written to the .dbf file. Alternatively a struct can be supplied with 43## attibute names contained in field "FieldName" (preferrably camelcased as 44## shown, but Octave treats this field's name as case-insensitive). If 45## argument @var{atts} is omitted all attributes will be written to the 46## shapefile. 47## 48## shapewrite produces 3 files, i.e. a .shp file (the actual shape file), 49## a .shx file (index file), and a .dbf file (dBase type 3) with the contents 50## of additional attribute fields other than Geometry, X/Y or Lat/Lon, M, Z, 51## and BoundingBox. If no additional attributes are present, a .dbf file is 52## created with the shape feature numbers as contents of field "ID". 53## 54## Return argument @var{status} is set to 1 if the complete shape file set was 55## written successfully, to 0 otherwise. 56## 57## Ref: http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf 58## 59## @seealso{shapedraw, shapeinfo, shaperead} 60## @end deftypefn 61 62## Author: Philip Nienhuis <prnienhuis@users.sf.net> 63## Created: 2014-12-30 64 65function [status] = shapewrite (shp, fname="", atts=[]) 66 67 persistent n_dbfspec = 0; 68 status = 0; 69 70 ## Input validation 71 if (nargin < 1) 72 print_usage; 73 endif 74 75 ## Assess shape variable type (oct or ml/geo ml/map) 76 if (! isstruct (shp)) 77 error ("shapewrite: [map-, geo-] struct expected for argument #1"); 78 else 79 ## Yep. Find out what type 80 fldn = fieldnames (shp); 81 if (ismember ("vals", fldn) && ismember ("shpbox", fldn)) 82 ## Assume it is an Octave-style struct read by shaperead 83 otype = 0; 84 warning ("shapewrite: only Matlab-type map/geostructs can be written\n"); 85 return; 86 elseif (ismember ("Geometry", fldn) && all (ismember ({"X", "Y"}, fldn))) 87 ## Assume it is a Matlab-style mapstruct 88 otype = 1; 89 elseif (ismember ("Geometry", fldn) && all (ismember ({"Lat", "Lon"}, fldn))) 90 ## Assume it is a Matlab-style geostruct 91 otype = 2; 92 else 93 ## Not a supported struct type 94 error ("shapewrite: unsupported struct type.\n") 95 endif 96 endif 97 98 ## FIXME struct field type validation 99 100 if (strcmpi (shp(1).Geometry, "MultiPatch")) 101 error ("shapewrite: MultiPatch type not supported"); 102 endif 103 104 ## Check for dbfwrite function 105 if (isempty (which ("dbfwrite"))) 106 error ("shapewrite: dbfwrite function not found. (io package installed \ 107and loaded?)"); 108 return; 109 endif 110 111 ## Check file name 112 if (isempty (fname)) 113 error ("shapewrite: filename expected for input argument #2"); 114 else 115 [pth, fnm, ext] = fileparts (fname); 116 if (isempty (ext)) 117 bname = fname; 118 fname = [fname ".shp"]; 119 ## Later on bname.shx and bname.dbf will be created 120 elseif (isempty (pth)) 121 bname = fnm; 122 else 123 bname = [pth filesep fnm]; 124 endif 125 endif 126 127 ## Check optional 3rd argument 128 if (nargin > 2) 129 if (isstruct (atts)) 130 if (! n_dbfspec) 131 warning ("shapewrite: DbfSpec not implemented; including requested \ 132attributes\n"); 133 n_dbfspec = 1; 134 endif 135 ## Get attribute names from field "FieldName"; allow lowercase and camelcase 136 ## Index of case-insensitive matches of "FieldName" 137 fnidx = find (ismember ("fieldname", lower (fieldnames (atts)))); 138 if (! isempty (fnidx)) 139 ## Get fieldnames out of first match 140 atts = {atts.(fieldnames (atts (fnidx)){1})}; 141 else 142 warning ("shapewrite: no field 'fieldname' (case-insensitive) found \ 143in struct\n=> input arg. #3 ignored"); 144 endif 145 elseif (! iscellstr (atts) && ! ischar (atts)) 146 error ("shapewrite: arg.#3: attribute name or cellstr array of attribute names expected"); 147 endif 148 ## Check if requested attributes exist at all in shapestruct 149 atts = unique (atts); 150 mtch = ! ismember (atts, fieldnames (shp)); 151 if (any (mtch)) 152 warning ("shapewrite: requested attribute(s) '%s' not in shapestruct\n", ... 153 strjoin (atts(find (mtch)), "', '")); 154 atts(mtch) = []; 155 endif 156 endif 157 158 ## Prepare a few things 159 numfeat = numel (shp); 160 if (abs (otype) >= 1) 161 stype = find (ismember ({"point", "multipoint", "line", "polygon"}, ... 162 lower (shp(1).Geometry))); 163## strmatch (lower (shp(1).Geometry), ... 164## {"point", "multipoint", "line", "polygon"}); 165 if (isempty (stype)) 166 ## Not a supported struct type 167 error ("shapewrite: unsupported struct type.\n") 168 else 169 stype = [1, 8, 3, 5](stype); 170 endif 171 172 ## Preprocess geostructs 173 if (abs (otype) == 2) 174 ## Change Lat/Lon fields into X/Y 175 [shp.X] = deal (shp.Lon); 176 [shp.Y] = deal (shp.Lat); 177 endif 178 ## "Point" need not have a BoundingBox field => add a dummy if not found 179 if (stype == 1 && ! ismember ("BoundingBox", fldn)) 180 [shp.BoundingBox] = deal ([0, 0; 0, 0]); 181 endif 182 if (any (ismember ({"M", "Z"}, fldn))) 183 if (ismember ({"Z"}, fldn)) 184 stype += 10; 185 else 186 stype += 20; 187 endif 188 otype = -otype; 189 endif 190 endif 191 192 ## Only now (after input checks) open .shp and .shx files & rewind just to be sure 193 fids = fopen (fname, "w"); 194 if (fids < 0) 195 error ("shapewrite: shapefile %s can't be opened for writing\n", fname); 196 endif 197 fseek (fids, 0, "bof"); 198 fidx = fopen ([ bname ".shx" ], "w"); 199 if (fidx < 0) 200 error ("shapewrite: index file %s can't be opened for writing\n", fname); 201 endif 202 fseek (fidx, 0, "bof"); 203 204 ## Write headers in .shp & .shx (identical). First magic number 9994 + six 205 ## zeros, the last zero in .shp is a placeholder for the yet unknown .shp 206 ## file length. 207 fwrite (fids, [9994 0 0 0 0 0 0], "int32", 0, "ieee-be"); 208 fwrite (fidx, [9994 0 0 0 0 0], "int32", 0, "ieee-be"); 209 ## For .shx the file length in 16-bit words (single) is known: 210 fwrite (fidx, ((numfeat * 4) + 50), "int32", 0, "ieee-be"); 211 ## Next, shp file version 212 fwrite (fids, 1000, "int32"); 213 fwrite (fidx, 1000, "int32"); 214 ## Shape feature type 215 fwrite (fids, stype, "int32"); 216 fwrite (fidx, stype, "int32"); 217 ## Bounding box. Can be run later for ML type shape structs. Fill with zeros 218 fwrite (fids, [0 0 0 0 0 0 0 0], "double"); 219 fwrite (fidx, [0 0 0 0 0 0 0 0], "double"); 220 ## Prepare BoundingBox limits 221 xMin = yMin = zMin = mMin = Inf; 222 xMax = yMax = zMax = mMax = -Inf; 223 224 ## Skip to start of first record position 225% fseek (fids, 100, "bof"); 226% fseek (fidx, 100, "bof"); 227 228 ## Write shape features one by one 229 if (abs (otype) >= 1) 230 for ishp=1:numfeat 231 ## Write record start pos to .shx file 232 fwrite (fidx, ftell (fids) / 2, "int32", 0, "ieee-be"); 233 234 ## Prepare multipart polygons/lines. 235 ## Find pointers to separators 236 idx = [ 0 find(isnan (shp(ishp).X)) ]; 237 ## Eliminate trailing NaN rows 238 if (isnan (shp(ishp).X)) 239 idx(end) = []; 240 endif 241 ## Augment idx for later on 242 idx = unique ([ idx (numel (shp(ishp).X)+1) ]); 243 ## Remove NaN separators 244 idn = find (! isfinite (shp(ishp).X)); 245 shp(ishp).X(idn) = []; 246 shp(ishp).Y(idn) = []; 247 if (stype >= 10) 248 shp(ishp).Z(idn) = []; 249 endif 250 if (stype >= 10) 251 shp(ishp).M(idn) = []; 252 ## M as need not be present (we assume a NaN value then). 253 ## Set M to a value less than -1e39 254 idm = find (! isfinite (shp(ishp).M)); 255 shp(ishp).M(idm) = -1.1e39; 256 endif 257 ## Nr. of vertices 258 npt = numel (shp(ishp).X); 259 ## Pointers to parts 260 nptr = idx(1:end-1) .- [0:numel(idx)-2]; 261 262 ## Write record contents 263 switch (stype) 264 case 1 ## Point 265 reclen = 10; 266 ## Write record index number & content length (fixed) 267 fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); 268 fwrite (fidx, reclen, "int32", 0, "ieee-be"); 269 ## Write shape type 270 fwrite (fids, stype, "int32"); 271 ## Simply write XY cordinates 272 fwrite (fids, [shp(ishp).X shp(ishp).Y], "double"); 273 ## Update overall BoundingBox 274 xMin = min (xMin, shp(ishp).X); 275 xMax = max (xMax, shp(ishp).X); 276 yMin = min (yMin, shp(ishp).Y); 277 yMax = max (yMax, shp(ishp).Y); 278 279 case 11 ## PointZ 280 reclen = 18; 281 ## Write record index number & content length (fixed) 282 fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); 283 fwrite (fidx, reclen, "int32", 0, "ieee-be"); 284 ## Write shape type 285 fwrite (fids, stype, "int32"); 286 ## Simply write XY coordinates & Z & M value 287 fwrite (fids, [shp(ishp).X shp(ishp).Y shp(ishp).Z shp(ishp).M], "double"); 288 ## Update overall BoundingBox 289 xMin = min (xMin, shp(ishp).X); 290 xMax = max (xMax, shp(ishp).X); 291 yMin = min (yMin, shp(ishp).Y); 292 yMax = max (yMax, shp(ishp).Y); 293 zMin = min (zMin, shp(ishp).Z); 294 zMax = max (zMax, shp(ishp).Z); 295 mMin = min (mMin, shp(ishp).M); 296 mMax = max (mMax, shp(ishp).M); 297 298 case 21 ## PointM 299 reclen = 14; 300 ## Write record index number & content length (fixed) 301 fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); 302 fwrite (fidx, reclen, "int32", 0, "ieee-be"); 303 ## Write shape type 304 fwrite (fids, stype, "int32"); 305 ## Simply write XY coordinates & M-value 306 fwrite (fids, [shp(ishp).X shp(ishp).Y shp(ishp).M], "double"); 307 ## Update overall BoundingBox 308 xMin = min (xMin, shp(ishp).X); 309 xMax = max (xMax, shp(ishp).X); 310 yMin = min (yMin, shp(ishp).Y); 311 yMax = max (yMax, shp(ishp).Y); 312 mMin = min (mMin, shp(ishp).M); 313 mMax = max (mMax, shp(ishp).M); 314 315 case 8 ## MultiPoint 316 reclen = (40 + 16 * numel (shp(ishp).X)) / 2; 317 ## Write record index number & content length (fixed) 318 fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); 319 fwrite (fidx, reclen, "int32", 0, "ieee-be"); 320 ## Write shape type (+4) 321 fwrite (fids, stype, "int32"); 322 ## Write bounding box (+32 -> 36) 323 fwrite (fids, [shp(ishp).BoundingBox'(:)]', "double"); 324 ## Update overall BoundingBox 325 xMin = min (xMin, min (shp(ishp).X)); 326 xMax = max (xMax, max (shp(ishp).X)); 327 yMin = min (yMin, min (shp(ishp).Y)); 328 yMax = max (yMax, max (shp(ishp).Y)); 329 ## Write nr of points and XY data (+4 -> 40 + Nx16) 330 fwrite (fids, numel (shp(ishp).X), "int32"); 331 fwrite (fids, [shp(ishp).X' shp(ishp).Y']', "double"); 332 333 case 18 ## MultiPointZ 334 reclen = (72 + 32 * numel (shp(ishp).X)) / 2; 335 ## Write record index number & content length (fixed) 336 fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); 337 fwrite (fidx, reclen, "int32", 0, "ieee-be"); 338 ## Write shape type 339 fwrite (fids, stype, "int32"); 340 ## Write bounding box 341 fwrite (fids, [shp(ishp).BoundingBox'(:)]', "double"); 342 ## Update overall BoundingBox 343 xMin = min (xMin, min (shp(ishp).X)); 344 xMax = max (xMax, max (shp(ishp).X)); 345 yMin = min (yMin, min (shp(ishp).Y)); 346 yMax = max (yMax, max (shp(ishp).Y)); 347 zMin = min (zMin, min (shp(ishp).Z)); 348 zMax = max (zMax, max (shp(ishp).Z)); 349 mMin = min (mMin, min (shp(ishp).M)); 350 mMax = max (mMax, max (shp(ishp).M)); 351 ## Write Nr of points and XY data 352 fwrite (fids, numel (shp(ishp).X), "int32"); 353 fwrite (fids, [shp(ishp).X' shp(ishp).Y']', "double"); 354 ## Write Z/M range and -data in turn 355 fwrite (fids, [min(shp(ishp).Z) max(shp(ishp).Z) shp(ishp).Z], "double"); 356 fwrite (fids, [min(shp(ishp).M) max(shp(ishp).M) shp(ishp).M], "double"); 357 358 case 28 ## MultiPointM 359 reclen = (56 + 24 * numel (shp(ishp).X)) / 2; 360 ## Write record index number & content length (fixed) 361 fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); 362 fwrite (fidx, reclen, "int32", 0, "ieee-be"); 363 ## Write shape type 364 fwrite (fids, stype, "int32"); 365 ## Write bounding box 366 fwrite (fids, [shp(ishp).BoundingBox'(:)]', "double"); 367 ## Update overall BoundingBox 368 xMin = min (xMin, min (shp(ishp).X)); 369 xMax = max (xMax, max (shp(ishp).X)); 370 yMin = min (yMin, min (shp(ishp).Y)); 371 yMax = max (yMax, max (shp(ishp).Y)); 372 mMin = min (mMin, min (shp(ishp).M)); 373 mMax = max (mMax, max (shp(ishp).M)); 374 ## Write Nr of points and XY data 375 fwrite (fids, numel (shp(ishp).X), "int32"); 376 fwrite (fids, [shp(ishp).X' shp(ishp).Y']', "double"); 377 ## Write M range and M data 378 fwrite (fids, [min(shp(ishp).M) max(shp(ishp).M) shp(ishp).M], "double"); 379 380 case {3, 5} ## Polyline/-gon 381 ## Content length 382 reclen = (44 + (numel(idx)-1)*4 + 2*8*npt) / 2; 383 ## Write record index number & content length (fixed) 384 fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); 385 fwrite (fidx, reclen, "int32", 0, "ieee-be"); 386 ## Write shape type 387 fwrite (fids, stype, "int32"); 388 ## Write bounding box 389 fwrite (fids, [shp(ishp).BoundingBox'(:)]', "double"); 390 ## Update overall BoundingBox 391 xMin = min (xMin, min (shp(ishp).X)); 392 xMax = max (xMax, max (shp(ishp).X)); 393 yMin = min (yMin, min (shp(ishp).Y)); 394 yMax = max (yMax, max (shp(ishp).Y)); 395 ## Write number of parts, number of points, part pointers 396 fwrite (fids, [(numel(idx)-1) npt nptr ], "int32"); 397 fwrite (fids, [shp(ishp).X' shp(ishp).Y']'(:), "double"); 398 399 case {13, 15} ## Polyline/-gonZ 400 ## Content length 401 reclen = (76 + (numel(idx)-1)*4 + 4*8*npt) / 2; 402 ## Write record index number & content length (fixed) 403 fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); 404 fwrite (fidx, reclen, "int32", 0, "ieee-be"); 405 ## Shape type 406 fwrite (fids, stype, "int32"); 407 ## Write bounding box 408 fwrite (fids, [shp(ishp).BoundingBox'(:)]', "double"); 409 ## Update overall BoundingBox 410 xMin = min (xMin, min (shp(ishp).X)); 411 xMax = max (xMax, max (shp(ishp).X)); 412 yMin = min (yMin, min (shp(ishp).Y)); 413 yMax = max (yMax, max (shp(ishp).Y)); 414 zMin = min (zMin, min (shp(ishp).Z)); 415 zMax = max (zMax, max (shp(ishp).Z)); 416 mMin = min (mMin, min (shp(ishp).M)); 417 mMax = max (mMax, max (shp(ishp).M)); 418 ## Write number of parts, number of points, part pointers 419 fwrite (fids, [(numel(idx)-1) npt nptr ], "int32"); 420 ## Write XY data 421 fwrite (fids, [shp(ishp).X' shp(ishp).Y']'(:), "double"); 422 fwrite (fids, [min(shp(ishp).Z) max(shp(ishp).Z) ... 423 shp(ishp).Z], "double"); 424 fwrite (fids, [min(shp(ishp).M) max(shp(ishp).M) ... 425 shp(ishp).M], "double"); 426 427 case {23, 25} ## Polyline/-gonM 428 ## Content length 429 reclen = (60 + (numel(idx)-1)*4 + 3*8*npt) / 2; 430 ## Write record index number & content length (fixed) 431 fwrite (fids, [ishp reclen], "int32", 0, "ieee-be"); 432 fwrite (fidx, reclen, "int32", 0, "ieee-be"); 433 ## Write shape type 434 fwrite (fids, stype, "int32"); 435 ## Write bounding box 436 fwrite (fids, [shp(ishp).BoundingBox'(:)]', "double"); 437 ## Update overall BoundingBox 438 xMin = min (xMin, min (shp(ishp).X)); 439 xMax = max (xMax, max (shp(ishp).X)); 440 yMin = min (yMin, min (shp(ishp).Y)); 441 yMax = max (yMax, max (shp(ishp).Y)); 442 mMin = min (mMin, min (shp(ishp).M)); 443 mMax = max (mMax, max (shp(ishp).M)); 444 ## Write number of parts, number of points, part pointers 445 fwrite (fids, [(numel(idx)-1) npt nptr ], "int32"); 446 ## Write XY data 447 fwrite (fids, [shp(ishp).X' shp(ishp).Y']'(:), "double"); 448 ## Write M range and M data 449 fwrite (fids, [min(shp(ishp).M) max(shp(ishp).M) ... 450 shp(ishp).M], "double"); 451 452 otherwise 453 ## Future shape types or types unsupported yet (MultiPatch) 454 455 endswitch 456 endfor 457 endif 458 459 ## Write file length and overall BoundingBox into .shp header 460 flen = ftell (fids); 461 fseek (fids, 24, "bof"); 462 fwrite (fids, flen/2, "int32", 0, "ieee-be"); 463 fseek (fids, 36, "bof"); 464 fwrite (fids, [xMin yMin xMax yMax], "double"); 465 ## Same for .shx header 466 xlen = ftell (fidx); 467 fseek (fidx, 24, "bof"); 468 fwrite (fidx, xlen/2, "int32", 0, "ieee-be"); 469 fseek (fidx, 36, "bof"); 470 fwrite (fidx, [xMin yMin xMax yMax], "double"); 471 if (stype > 10) 472 ## +-Inf & NaN not allowed in shapefiles 473 zm = [zMin zMax mMin mMax]; 474 zm (! isfinite (zm)) = -1e-39; 475 fwrite (fids, zm, "double"); 476 fwrite (fidx, zm, "double"); 477 endif 478 479 ## Close files 480 fclose (fids); 481 fclose (fidx); 482 483 ## Write .dbf file. 484 ## Remove basic attributes 485 if (abs (otype) == 1) 486 ## Attributes + shp data in mapstruct 487 shp = rmfield (shp, {"Geometry", "BoundingBox", "X", "Y"}); 488 elseif (abs (otype) == 2) 489 ## Attributes + shp data in geostruct 490 shp = rmfield (shp, {"Geometry", "BoundingBox", "Lat", "Lon", "X", "Y"}); 491 endif 492 if (otype < 1) 493 shp = rmfield (shp, {"M"}); 494 if (isfield (shp, "Z")) 495 shp = rmfield (shp, {"Z"}); 496 endif 497 endif 498 499 ## Write rest of attributes 500 if (nargin == 3) 501 ## Only write user-specified attribute selection 502 fldn = fieldnames (shp); 503 ## First remove regular attributes (with a value for each vertex) 504 fldn (ismember (fldn, {"Geometry", "BoundingBox", "Lat", "Lon", "X", "Y", "M", "Z"})) = []; 505 ## Next, only retain user-specified attributes 506 shp = rmfield (shp, setdiff (fldn, atts)); 507 endif 508 attribs = cell (numfeat + 1, numel (fieldnames (shp))); 509 if (! isempty (attribs)) 510 attribs(1, :) = fieldnames (shp); 511 attribs(2:end, :) = (squeeze (struct2cell (shp)))'; 512 else 513 ## Substitute ID attribute 514 attribs{1, 1} = "ID"; 515 [attribs{2:end}] = deal (num2cell ([1:size(shp, 2)]){:}); 516 endif 517 try 518 status = dbfwrite ([ bname ".dbf"], attribs); 519 status = 1; 520 catch 521 warning ("shapewrite: writing attributes to file %s failed\n", [bname ".dbf"]); 522 end_try_catch 523 524endfunction 525 526 527## Test various geometries: (1) Point 528%!test 529%! shp.Geometry = "Point"; 530%! shp.X = 10; 531%! shp.Y = 20; 532%! shp.Z = 30; 533%! shp.M = -1; 534%! shp.attr_1 = "Attribute1"; 535%! shp.attr_Z = "AttributeA"; 536%! shp(2).Geometry = "Point"; 537%! shp(2).X = 11; 538%! shp(2).Y = 25; 539%! shp(2).Z = 32; 540%! shp(2).M = -2; 541%! shp(2).attr_1 = "Attribute2"; 542%! shp(2).attr_Z = "AttributeB"; 543%! fname = tempname (); 544%! sts = shapewrite (shp, fname); 545%! assert (sts, 1, eps); 546%! ## Check index file 547%%! fx = fopen ([fname ".shx"], "r"); 548%%! fseek (fx, 100, "bof"); 549%%! shxinfo = fread (fx, "*int32", 0, "ieee-be"); 550%%! assert (shxinfo, int32 ([50; 36; 72; 36])); 551%%! fclose (fx); 552%! ## Check on filesizes, based on Esri shapewrite doc 553%! assert (stat ([fname ".shp"]).size, 188, 1e-10); 554%! assert (stat ([fname ".shx"]).size, 116, 1e-10); 555%! shp2 = shaperead ([fname ".shp"]); 556%! assert (size (shp2), [2 1]); 557%! flds = fieldnames (shp2); 558%! fields = {"Geometry", "X", "Y", "attr_1", "attr_Z"}; 559%! ism = ismember (fields, flds); 560%! ## Do we have only those fields? 561%! assert (numel (ism), numel (fields)); 562%! ## Do we have only those fields? 563%! assert (sum (ism), numel (ism)); 564%! unlink ([fname ".shp"]); 565%! unlink ([fname ".shx"]); 566%! assert ([shp2.X shp2.Y], [10 11 20 25], 1e-10); 567%! assert ({shp2.Geometry}, {"Point", "Point"}); 568%! assert ({shp2.attr_1}, {"Attribute1", "Attribute2"}); 569 570## Test various geometries: (2) Line & Polygon 571%!test 572%! shp.Geometry = "Line"; 573%! shp.BoundingBox = [9 110; 19 120]; 574%! shp.X = [10 110 NaN 9 109]; 575%! shp.Y = [20 120 NaN 19 119]; 576%! shp.Z = [30 130 NaN 29 129]; 577%! shp.M = [-1 1 NaN -11 11]; 578%! shp.attr_1 = "Attribute1"; 579%! shp.attr_Z = "AttributeA"; 580%! shp(2).Geometry = "Line"; 581%! shp(2).BoundingBox = [11 211; 24 225]; 582%! shp(2).X = [11 111 211 NaN 11 110 NaN 12 200]; 583%! shp(2).Y = [25 125 225 NaN 24 124 NaN 26 200]; 584%! shp(2).Z = [32 132 232 NaN 31 131 NaN 33 200]; 585%! shp(2).M = [-2 NaN -3 NaN -22 22 NaN -33 33]; 586%! shp(2).attr_1 = "Attribute2"; 587%! shp(2).attr_Z = "AttributeB"; 588%! fname = tempname (); 589%! sts = shapewrite (shp, fname); 590%! assert (sts, 1, eps); 591%! ## Check index file 592%%! fx = fopen ([fname ".shx"], "r"); 593%%! fseek (fx, 100, "bof"); 594%%! shxinfo = fread (fx, "*int32", 0, "ieee-be"); 595%%! assert (shxinfo, int32 ([50; 106; 160; 156])); 596%%! fclose (fx); 597%! ## Check on filesizes, based on Esri shapewrite doc 598%! assert (stat ([fname ".shp"]).size, 640, 1e-10); 599%! assert (stat ([fname ".shx"]).size, 116, 1e-10); 600%! 601%! ## Check Matlab-style reading 602%! shp2 = shaperead ([fname ".shp"]); 603%! assert (size (shp2), [2 1]); 604%! flds = fieldnames (shp2); 605%! fields = {"Geometry", "BoundingBox", "X", "Y", "attr_1", "attr_Z"}; 606%! ism = ismember (fields, flds); 607%! ## Do we have only those fields? 608%! assert (numel (ism), numel (fields)); 609%! ## Do we have only those fields? 610%! assert (sum (ism), numel (ism)); 611%! assert ([shp2.BoundingBox], [9 110 11 211; 19 120 24 225], 1e-10); 612%! assert ({shp2.attr_Z}, {"AttributeA", "AttributeB"}); 613%! 614%! ## Re-read using M & Z values, read only 2nd feature 615%! shp2 = shaperead ([fname ".shp"], "e", "rec", 2); 616%! assert (size (shp2), [1 1]); 617%! flds = fieldnames (shp2); 618%! fields = {"Geometry", "BoundingBox", "X", "Y", "Z", "M", "attr_1", "attr_Z"}; 619%! ism = ismember (fields, flds); 620%! ## Do we have only those fields? 621%! assert (numel (ism), numel (fields)); 622%! ## Do we have only those fields? 623%! assert (sum (ism), numel (ism)); 624%! ## Check regular numerical values 625%! assert ([shp2.X; shp2.Y; shp2.Z], ... 626%! [11, 111, 211, NaN, 11, 110, NaN, 12, 200; ... 627%! 25, 125, 225, NaN, 24, 124, NaN, 26, 200; ... 628%! 32, 132, 232, NaN, 31, 131, NaN, 33, 200], 1e-10); 629%! ## Check missing M value 630%! assert (shp2.M, [-2, NaN, -3, NaN, -22, 22, NaN, -33, 33], 1e-10); 631%! assert ({shp2.attr_1, shp2.attr_Z}, {"Attribute2", "AttributeB"}); 632%! 633%! unlink ([fname ".shp"]); 634%! unlink ([fname ".shx"]); 635