1""" 2shapefile.py 3Provides read and write support for ESRI Shapefiles. 4author: jlawhead<at>geospatialpython.com 5version: 2.1.3 6Compatible with Python versions 2.7-3.x 7""" 8 9__version__ = "2.1.3" 10 11from struct import pack, unpack, calcsize, error, Struct 12import os 13import sys 14import time 15import array 16import tempfile 17import warnings 18import io 19from datetime import date 20 21 22# Constants for shape types 23NULL = 0 24POINT = 1 25POLYLINE = 3 26POLYGON = 5 27MULTIPOINT = 8 28POINTZ = 11 29POLYLINEZ = 13 30POLYGONZ = 15 31MULTIPOINTZ = 18 32POINTM = 21 33POLYLINEM = 23 34POLYGONM = 25 35MULTIPOINTM = 28 36MULTIPATCH = 31 37 38SHAPETYPE_LOOKUP = { 39 0: 'NULL', 40 1: 'POINT', 41 3: 'POLYLINE', 42 5: 'POLYGON', 43 8: 'MULTIPOINT', 44 11: 'POINTZ', 45 13: 'POLYLINEZ', 46 15: 'POLYGONZ', 47 18: 'MULTIPOINTZ', 48 21: 'POINTM', 49 23: 'POLYLINEM', 50 25: 'POLYGONM', 51 28: 'MULTIPOINTM', 52 31: 'MULTIPATCH'} 53 54TRIANGLE_STRIP = 0 55TRIANGLE_FAN = 1 56OUTER_RING = 2 57INNER_RING = 3 58FIRST_RING = 4 59RING = 5 60 61PARTTYPE_LOOKUP = { 62 0: 'TRIANGLE_STRIP', 63 1: 'TRIANGLE_FAN', 64 2: 'OUTER_RING', 65 3: 'INNER_RING', 66 4: 'FIRST_RING', 67 5: 'RING'} 68 69 70# Python 2-3 handling 71 72PYTHON3 = sys.version_info[0] == 3 73 74if PYTHON3: 75 xrange = range 76 izip = zip 77else: 78 from itertools import izip 79 80 81# Helpers 82 83MISSING = [None,''] 84NODATA = -10e38 # as per the ESRI shapefile spec, only used for m-values. 85 86if PYTHON3: 87 def b(v, encoding='utf-8', encodingErrors='strict'): 88 if isinstance(v, str): 89 # For python 3 encode str to bytes. 90 return v.encode(encoding, encodingErrors) 91 elif isinstance(v, bytes): 92 # Already bytes. 93 return v 94 elif v is None: 95 # Since we're dealing with text, interpret None as "" 96 return b"" 97 else: 98 # Force string representation. 99 return str(v).encode(encoding, encodingErrors) 100 101 def u(v, encoding='utf-8', encodingErrors='strict'): 102 if isinstance(v, bytes): 103 # For python 3 decode bytes to str. 104 return v.decode(encoding, encodingErrors) 105 elif isinstance(v, str): 106 # Already str. 107 return v 108 elif v is None: 109 # Since we're dealing with text, interpret None as "" 110 return "" 111 else: 112 # Force string representation. 113 return bytes(v).decode(encoding, encodingErrors) 114 115 def is_string(v): 116 return isinstance(v, str) 117 118else: 119 def b(v, encoding='utf-8', encodingErrors='strict'): 120 if isinstance(v, unicode): 121 # For python 2 encode unicode to bytes. 122 return v.encode(encoding, encodingErrors) 123 elif isinstance(v, bytes): 124 # Already bytes. 125 return v 126 elif v is None: 127 # Since we're dealing with text, interpret None as "" 128 return "" 129 else: 130 # Force string representation. 131 return unicode(v).encode(encoding, encodingErrors) 132 133 def u(v, encoding='utf-8', encodingErrors='strict'): 134 if isinstance(v, bytes): 135 # For python 2 decode bytes to unicode. 136 return v.decode(encoding, encodingErrors) 137 elif isinstance(v, unicode): 138 # Already unicode. 139 return v 140 elif v is None: 141 # Since we're dealing with text, interpret None as "" 142 return u"" 143 else: 144 # Force string representation. 145 return bytes(v).decode(encoding, encodingErrors) 146 147 def is_string(v): 148 return isinstance(v, basestring) 149 150 151# Begin 152 153class _Array(array.array): 154 """Converts python tuples to lists of the appropritate type. 155 Used to unpack different shapefile header parts.""" 156 def __repr__(self): 157 return str(self.tolist()) 158 159def signed_area(coords): 160 """Return the signed area enclosed by a ring using the linear time 161 algorithm. A value >= 0 indicates a counter-clockwise oriented ring. 162 """ 163 xs, ys = map(list, list(zip(*coords))[:2]) # ignore any z or m values 164 xs.append(xs[1]) 165 ys.append(ys[1]) 166 return sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(coords)))/2.0 167 168def ring_bbox(coords): 169 """Calculates and returns the bounding box of a ring. 170 """ 171 xs,ys = zip(*coords) 172 bbox = min(xs),min(ys),max(xs),max(ys) 173 return bbox 174 175def bbox_overlap(bbox1, bbox2): 176 """Tests whether two bounding boxes overlap, returning a boolean 177 """ 178 xmin1,ymin1,xmax1,ymax1 = bbox1 179 xmin2,ymin2,xmax2,ymax2 = bbox2 180 overlap = (xmin1 <= xmax2 and xmax1 >= xmin2 and ymin1 <= ymax2 and ymax1 >= ymin2) 181 return overlap 182 183def bbox_contains(bbox1, bbox2): 184 """Tests whether bbox1 fully contains bbox2, returning a boolean 185 """ 186 xmin1,ymin1,xmax1,ymax1 = bbox1 187 xmin2,ymin2,xmax2,ymax2 = bbox2 188 contains = (xmin1 < xmin2 and xmax1 > xmax2 and ymin1 < ymin2 and ymax1 > ymax2) 189 return contains 190 191def ring_contains_point(coords, p): 192 """Fast point-in-polygon crossings algorithm, MacMartin optimization. 193 194 Adapted from code by Eric Haynes 195 http://www.realtimerendering.com/resources/GraphicsGems//gemsiv/ptpoly_haines/ptinpoly.c 196 197 Original description: 198 Shoot a test ray along +X axis. The strategy, from MacMartin, is to 199 compare vertex Y values to the testing point's Y and quickly discard 200 edges which are entirely to one side of the test ray. 201 """ 202 tx,ty = p 203 204 # get initial test bit for above/below X axis 205 vtx0 = coords[0] 206 yflag0 = ( vtx0[1] >= ty ) 207 208 inside_flag = False 209 for vtx1 in coords[1:]: 210 yflag1 = ( vtx1[1] >= ty ) 211 # check if endpoints straddle (are on opposite sides) of X axis 212 # (i.e. the Y's differ); if so, +X ray could intersect this edge. 213 if yflag0 != yflag1: 214 xflag0 = ( vtx0[0] >= tx ) 215 # check if endpoints are on same side of the Y axis (i.e. X's 216 # are the same); if so, it's easy to test if edge hits or misses. 217 if xflag0 == ( vtx1[0] >= tx ): 218 # if edge's X values both right of the point, must hit 219 if xflag0: 220 inside_flag = not inside_flag 221 else: 222 # compute intersection of pgon segment with +X ray, note 223 # if >= point's X; if so, the ray hits it. 224 if ( vtx1[0] - (vtx1[1]-ty) * ( vtx0[0]-vtx1[0]) / (vtx0[1]-vtx1[1]) ) >= tx: 225 inside_flag = not inside_flag 226 227 # move to next pair of vertices, retaining info as possible 228 yflag0 = yflag1 229 vtx0 = vtx1 230 231 return inside_flag 232 233def ring_sample(coords, ccw=False): 234 """Return a sample point guaranteed to be within a ring, by efficiently 235 finding the first centroid of a coordinate triplet whose orientation 236 matches the orientation of the ring and passes the point-in-ring test. 237 The orientation of the ring is assumed to be clockwise, unless ccw 238 (counter-clockwise) is set to True. 239 """ 240 triplet = [] 241 def itercoords(): 242 # iterate full closed ring 243 for p in coords: 244 yield p 245 # finally, yield the second coordinate to the end to allow checking the last triplet 246 yield coords[1] 247 248 for p in itercoords(): 249 # add point to triplet (but not if duplicate) 250 if p not in triplet: 251 triplet.append(p) 252 253 # new triplet, try to get sample 254 if len(triplet) == 3: 255 # check that triplet does not form a straight line (not a triangle) 256 is_straight_line = (triplet[0][1] - triplet[1][1]) * (triplet[0][0] - triplet[2][0]) == (triplet[0][1] - triplet[2][1]) * (triplet[0][0] - triplet[1][0]) 257 if not is_straight_line: 258 # get triplet orientation 259 closed_triplet = triplet + [triplet[0]] 260 triplet_ccw = signed_area(closed_triplet) >= 0 261 # check that triplet has the same orientation as the ring (means triangle is inside the ring) 262 if ccw == triplet_ccw: 263 # get triplet centroid 264 xs,ys = zip(*triplet) 265 xmean,ymean = sum(xs) / 3.0, sum(ys) / 3.0 266 # check that triplet centroid is truly inside the ring 267 if ring_contains_point(coords, (xmean,ymean)): 268 return xmean,ymean 269 270 # failed to get sample point from this triplet 271 # remove oldest triplet coord to allow iterating to next triplet 272 triplet.pop(0) 273 274 else: 275 raise Exception('Unexpected error: Unable to find a ring sample point.') 276 277def ring_contains_ring(coords1, coords2): 278 '''Returns True if all vertexes in coords2 are fully inside coords1. 279 ''' 280 return all((ring_contains_point(coords1, p2) for p2 in coords2)) 281 282def organize_polygon_rings(rings): 283 '''Organize a list of coordinate rings into one or more polygons with holes. 284 Returns a list of polygons, where each polygon is composed of a single exterior 285 ring, and one or more interior holes. 286 287 Rings must be closed, and cannot intersect each other (non-self-intersecting polygon). 288 Rings are determined as exteriors if they run in clockwise direction, or interior 289 holes if they run in counter-clockwise direction. This method is used to construct 290 GeoJSON (multi)polygons from the shapefile polygon shape type, which does not 291 explicitly store the structure of the polygons beyond exterior/interior ring orientation. 292 ''' 293 # first iterate rings and classify as exterior or hole 294 exteriors = [] 295 holes = [] 296 for ring in rings: 297 # shapefile format defines a polygon as a sequence of rings 298 # where exterior rings are clockwise, and holes counterclockwise 299 if signed_area(ring) < 0: 300 # ring is exterior 301 exteriors.append(ring) 302 else: 303 # ring is a hole 304 holes.append(ring) 305 306 # if only one exterior, then all holes belong to that exterior 307 if len(exteriors) == 1: 308 # exit early 309 poly = [exteriors[0]] + holes 310 polys = [poly] 311 return polys 312 313 # multiple exteriors, ie multi-polygon, have to group holes with correct exterior 314 # shapefile format does not specify which holes belong to which exteriors 315 # so have to do efficient multi-stage checking of hole-to-exterior containment 316 elif len(exteriors) > 1: 317 # exit early if no holes 318 if not holes: 319 polys = [] 320 for ext in exteriors: 321 poly = [ext] 322 polys.append(poly) 323 return polys 324 325 # first determine each hole's candidate exteriors based on simple bbox contains test 326 hole_exteriors = dict([(hole_i,[]) for hole_i in xrange(len(holes))]) 327 exterior_bboxes = [ring_bbox(ring) for ring in exteriors] 328 for hole_i in hole_exteriors.keys(): 329 hole_bbox = ring_bbox(holes[hole_i]) 330 for ext_i,ext_bbox in enumerate(exterior_bboxes): 331 if bbox_contains(ext_bbox, hole_bbox): 332 hole_exteriors[hole_i].append( ext_i ) 333 334 # then, for holes with still more than one possible exterior, do more detailed hole-in-ring test 335 for hole_i,exterior_candidates in hole_exteriors.items(): 336 337 if len(exterior_candidates) > 1: 338 # get hole sample point 339 hole_sample = ring_sample(holes[hole_i], ccw=True) 340 # collect new exterior candidates 341 new_exterior_candidates = [] 342 for ext_i in exterior_candidates: 343 # check that hole sample point is inside exterior 344 hole_in_exterior = ring_contains_point(exteriors[ext_i], hole_sample) 345 if hole_in_exterior: 346 new_exterior_candidates.append(ext_i) 347 348 # set new exterior candidates 349 hole_exteriors[hole_i] = new_exterior_candidates 350 351 # if still holes with more than one possible exterior, means we have an exterior hole nested inside another exterior's hole 352 for hole_i,exterior_candidates in hole_exteriors.items(): 353 354 if len(exterior_candidates) > 1: 355 # exterior candidate with the smallest area is the hole's most immediate parent 356 ext_i = sorted(exterior_candidates, key=lambda x: abs(signed_area(exteriors[x])))[0] 357 hole_exteriors[hole_i] = [ext_i] 358 359 # separate out holes that are orphaned (not contained by any exterior) 360 orphan_holes = [] 361 for hole_i,exterior_candidates in list(hole_exteriors.items()): 362 if not exterior_candidates: 363 warnings.warn('Shapefile shape has invalid polygon: found orphan hole (not contained by any of the exteriors); interpreting as exterior.') 364 orphan_holes.append( hole_i ) 365 del hole_exteriors[hole_i] 366 continue 367 368 # each hole should now only belong to one exterior, group into exterior-holes polygons 369 polys = [] 370 for ext_i,ext in enumerate(exteriors): 371 poly = [ext] 372 # find relevant holes 373 poly_holes = [] 374 for hole_i,exterior_candidates in list(hole_exteriors.items()): 375 # hole is relevant if previously matched with this exterior 376 if exterior_candidates[0] == ext_i: 377 poly_holes.append( holes[hole_i] ) 378 poly += poly_holes 379 polys.append(poly) 380 381 # add orphan holes as exteriors 382 for hole_i in orphan_holes: 383 ext = holes[hole_i] # could potentially reverse their order, but in geojson winding order doesn't matter 384 poly = [ext] 385 polys.append(poly) 386 387 return polys 388 389 # no exteriors, be nice and assume due to incorrect winding order 390 else: 391 warnings.warn('Shapefile shape has invalid polygon: no exterior rings found (must have clockwise orientation); interpreting holes as exteriors.') 392 exteriors = holes # could potentially reverse their order, but in geojson winding order doesn't matter 393 polys = [[ext] for ext in exteriors] 394 return polys 395 396class Shape(object): 397 def __init__(self, shapeType=NULL, points=None, parts=None, partTypes=None): 398 """Stores the geometry of the different shape types 399 specified in the Shapefile spec. Shape types are 400 usually point, polyline, or polygons. Every shape type 401 except the "Null" type contains points at some level for 402 example vertices in a polygon. If a shape type has 403 multiple shapes containing points within a single 404 geometry record then those shapes are called parts. Parts 405 are designated by their starting index in geometry record's 406 list of shapes. For MultiPatch geometry, partTypes designates 407 the patch type of each of the parts. 408 """ 409 self.shapeType = shapeType 410 self.points = points or [] 411 self.parts = parts or [] 412 if partTypes: 413 self.partTypes = partTypes 414 415 @property 416 def __geo_interface__(self): 417 if self.shapeType in [POINT, POINTM, POINTZ]: 418 # point 419 if len(self.points) == 0: 420 # the shape has no coordinate information, i.e. is 'empty' 421 # the geojson spec does not define a proper null-geometry type 422 # however, it does allow geometry types with 'empty' coordinates to be interpreted as null-geometries 423 return {'type':'Point', 'coordinates':tuple()} 424 else: 425 return { 426 'type': 'Point', 427 'coordinates': tuple(self.points[0]) 428 } 429 elif self.shapeType in [MULTIPOINT, MULTIPOINTM, MULTIPOINTZ]: 430 if len(self.points) == 0: 431 # the shape has no coordinate information, i.e. is 'empty' 432 # the geojson spec does not define a proper null-geometry type 433 # however, it does allow geometry types with 'empty' coordinates to be interpreted as null-geometries 434 return {'type':'MultiPoint', 'coordinates':[]} 435 else: 436 # multipoint 437 return { 438 'type': 'MultiPoint', 439 'coordinates': [tuple(p) for p in self.points] 440 } 441 elif self.shapeType in [POLYLINE, POLYLINEM, POLYLINEZ]: 442 if len(self.parts) == 0: 443 # the shape has no coordinate information, i.e. is 'empty' 444 # the geojson spec does not define a proper null-geometry type 445 # however, it does allow geometry types with 'empty' coordinates to be interpreted as null-geometries 446 return {'type':'LineString', 'coordinates':[]} 447 elif len(self.parts) == 1: 448 # linestring 449 return { 450 'type': 'LineString', 451 'coordinates': [tuple(p) for p in self.points] 452 } 453 else: 454 # multilinestring 455 ps = None 456 coordinates = [] 457 for part in self.parts: 458 if ps == None: 459 ps = part 460 continue 461 else: 462 coordinates.append([tuple(p) for p in self.points[ps:part]]) 463 ps = part 464 else: 465 coordinates.append([tuple(p) for p in self.points[part:]]) 466 return { 467 'type': 'MultiLineString', 468 'coordinates': coordinates 469 } 470 elif self.shapeType in [POLYGON, POLYGONM, POLYGONZ]: 471 if len(self.parts) == 0: 472 # the shape has no coordinate information, i.e. is 'empty' 473 # the geojson spec does not define a proper null-geometry type 474 # however, it does allow geometry types with 'empty' coordinates to be interpreted as null-geometries 475 return {'type':'Polygon', 'coordinates':[]} 476 else: 477 # get all polygon rings 478 rings = [] 479 for i in xrange(len(self.parts)): 480 # get indexes of start and end points of the ring 481 start = self.parts[i] 482 try: 483 end = self.parts[i+1] 484 except IndexError: 485 end = len(self.points) 486 487 # extract the points that make up the ring 488 ring = [tuple(p) for p in self.points[start:end]] 489 rings.append(ring) 490 491 # organize rings into list of polygons, where each polygon is defined as list of rings. 492 # the first ring is the exterior and any remaining rings are holes (same as GeoJSON). 493 polys = organize_polygon_rings(rings) 494 495 # return as geojson 496 if len(polys) == 1: 497 return { 498 'type': 'Polygon', 499 'coordinates': polys[0] 500 } 501 else: 502 return { 503 'type': 'MultiPolygon', 504 'coordinates': polys 505 } 506 507 else: 508 raise Exception('Shape type "%s" cannot be represented as GeoJSON.' % SHAPETYPE_LOOKUP[self.shapeType]) 509 510 @staticmethod 511 def _from_geojson(geoj): 512 # create empty shape 513 shape = Shape() 514 # set shapeType 515 geojType = geoj["type"] if geoj else "Null" 516 if geojType == "Null": 517 shapeType = NULL 518 elif geojType == "Point": 519 shapeType = POINT 520 elif geojType == "LineString": 521 shapeType = POLYLINE 522 elif geojType == "Polygon": 523 shapeType = POLYGON 524 elif geojType == "MultiPoint": 525 shapeType = MULTIPOINT 526 elif geojType == "MultiLineString": 527 shapeType = POLYLINE 528 elif geojType == "MultiPolygon": 529 shapeType = POLYGON 530 else: 531 raise Exception("Cannot create Shape from GeoJSON type '%s'" % geojType) 532 shape.shapeType = shapeType 533 534 # set points and parts 535 if geojType == "Point": 536 shape.points = [ geoj["coordinates"] ] 537 shape.parts = [0] 538 elif geojType in ("MultiPoint","LineString"): 539 shape.points = geoj["coordinates"] 540 shape.parts = [0] 541 elif geojType in ("Polygon"): 542 points = [] 543 parts = [] 544 index = 0 545 for i,ext_or_hole in enumerate(geoj["coordinates"]): 546 if i == 0 and not signed_area(ext_or_hole) < 0: 547 # flip exterior direction 548 ext_or_hole = list(reversed(ext_or_hole)) 549 elif i > 0 and not signed_area(ext_or_hole) >= 0: 550 # flip hole direction 551 ext_or_hole = list(reversed(ext_or_hole)) 552 points.extend(ext_or_hole) 553 parts.append(index) 554 index += len(ext_or_hole) 555 shape.points = points 556 shape.parts = parts 557 elif geojType in ("MultiLineString"): 558 points = [] 559 parts = [] 560 index = 0 561 for linestring in geoj["coordinates"]: 562 points.extend(linestring) 563 parts.append(index) 564 index += len(linestring) 565 shape.points = points 566 shape.parts = parts 567 elif geojType in ("MultiPolygon"): 568 points = [] 569 parts = [] 570 index = 0 571 for polygon in geoj["coordinates"]: 572 for i,ext_or_hole in enumerate(polygon): 573 if i == 0 and not signed_area(ext_or_hole) < 0: 574 # flip exterior direction 575 ext_or_hole = list(reversed(ext_or_hole)) 576 elif i > 0 and not signed_area(ext_or_hole) >= 0: 577 # flip hole direction 578 ext_or_hole = list(reversed(ext_or_hole)) 579 points.extend(ext_or_hole) 580 parts.append(index) 581 index += len(ext_or_hole) 582 shape.points = points 583 shape.parts = parts 584 return shape 585 586 @property 587 def shapeTypeName(self): 588 return SHAPETYPE_LOOKUP[self.shapeType] 589 590class _Record(list): 591 """ 592 A class to hold a record. Subclasses list to ensure compatibility with 593 former work and to reuse all the optimizations of the builtin list. 594 In addition to the list interface, the values of the record 595 can also be retrieved using the field's name. For example if the dbf contains 596 a field ID at position 0, the ID can be retrieved with the position, the field name 597 as a key, or the field name as an attribute. 598 599 >>> # Create a Record with one field, normally the record is created by the Reader class 600 >>> r = _Record({'ID': 0}, [0]) 601 >>> print(r[0]) 602 >>> print(r['ID']) 603 >>> print(r.ID) 604 """ 605 606 def __init__(self, field_positions, values, oid=None): 607 """ 608 A Record should be created by the Reader class 609 610 :param field_positions: A dict mapping field names to field positions 611 :param values: A sequence of values 612 :param oid: The object id, an int (optional) 613 """ 614 self.__field_positions = field_positions 615 if oid is not None: 616 self.__oid = oid 617 else: 618 self.__oid = -1 619 list.__init__(self, values) 620 621 def __getattr__(self, item): 622 """ 623 __getattr__ is called if an attribute is used that does 624 not exist in the normal sense. For example r=Record(...), r.ID 625 calls r.__getattr__('ID'), but r.index(5) calls list.index(r, 5) 626 :param item: The field name, used as attribute 627 :return: Value of the field 628 :raises: AttributeError, if item is not a field of the shapefile 629 and IndexError, if the field exists but the field's 630 corresponding value in the Record does not exist 631 """ 632 try: 633 index = self.__field_positions[item] 634 return list.__getitem__(self, index) 635 except KeyError: 636 raise AttributeError('{} is not a field name'.format(item)) 637 except IndexError: 638 raise IndexError('{} found as a field but not enough values available.'.format(item)) 639 640 def __setattr__(self, key, value): 641 """ 642 Sets a value of a field attribute 643 :param key: The field name 644 :param value: the value of that field 645 :return: None 646 :raises: AttributeError, if key is not a field of the shapefile 647 """ 648 if key.startswith('_'): # Prevent infinite loop when setting mangled attribute 649 return list.__setattr__(self, key, value) 650 try: 651 index = self.__field_positions[key] 652 return list.__setitem__(self, index, value) 653 except KeyError: 654 raise AttributeError('{} is not a field name'.format(key)) 655 656 def __getitem__(self, item): 657 """ 658 Extends the normal list item access with 659 access using a fieldname 660 661 For example r['ID'], r[0] 662 :param item: Either the position of the value or the name of a field 663 :return: the value of the field 664 """ 665 try: 666 return list.__getitem__(self, item) 667 except TypeError: 668 try: 669 index = self.__field_positions[item] 670 except KeyError: 671 index = None 672 if index is not None: 673 return list.__getitem__(self, index) 674 else: 675 raise IndexError('"{}" is not a field name and not an int'.format(item)) 676 677 def __setitem__(self, key, value): 678 """ 679 Extends the normal list item access with 680 access using a fieldname 681 682 For example r['ID']=2, r[0]=2 683 :param key: Either the position of the value or the name of a field 684 :param value: the new value of the field 685 """ 686 try: 687 return list.__setitem__(self, key, value) 688 except TypeError: 689 index = self.__field_positions.get(key) 690 if index is not None: 691 return list.__setitem__(self, index, value) 692 else: 693 raise IndexError('{} is not a field name and not an int'.format(key)) 694 695 @property 696 def oid(self): 697 """The index position of the record in the original shapefile""" 698 return self.__oid 699 700 def as_dict(self, date_strings=False): 701 """ 702 Returns this Record as a dictionary using the field names as keys 703 :return: dict 704 """ 705 dct = dict((f, self[i]) for f, i in self.__field_positions.items()) 706 if date_strings: 707 for k,v in dct.items(): 708 if isinstance(v, date): 709 dct[k] = '{:04d}{:02d}{:02d}'.format(v.year, v.month, v.day) 710 return dct 711 712 def __repr__(self): 713 return 'Record #{}: {}'.format(self.__oid, list(self)) 714 715 def __dir__(self): 716 """ 717 Helps to show the field names in an interactive environment like IPython. 718 See: http://ipython.readthedocs.io/en/stable/config/integrating.html 719 720 :return: List of method names and fields 721 """ 722 default = list(dir(type(self))) # default list methods and attributes of this class 723 fnames = list(self.__field_positions.keys()) # plus field names (random order if Python version < 3.6) 724 return default + fnames 725 726class ShapeRecord(object): 727 """A ShapeRecord object containing a shape along with its attributes. 728 Provides the GeoJSON __geo_interface__ to return a Feature dictionary.""" 729 def __init__(self, shape=None, record=None): 730 self.shape = shape 731 self.record = record 732 733 @property 734 def __geo_interface__(self): 735 return {'type': 'Feature', 736 'properties': self.record.as_dict(date_strings=True), 737 'geometry': None if self.shape.shapeType == NULL else self.shape.__geo_interface__} 738 739class Shapes(list): 740 """A class to hold a list of Shape objects. Subclasses list to ensure compatibility with 741 former work and to reuse all the optimizations of the builtin list. 742 In addition to the list interface, this also provides the GeoJSON __geo_interface__ 743 to return a GeometryCollection dictionary.""" 744 745 def __repr__(self): 746 return 'Shapes: {}'.format(list(self)) 747 748 @property 749 def __geo_interface__(self): 750 # Note: currently this will fail if any of the shapes are null-geometries 751 # could be fixed by storing the shapefile shapeType upon init, returning geojson type with empty coords 752 return {'type': 'GeometryCollection', 753 'geometries': [shape.__geo_interface__ for shape in self]} 754 755class ShapeRecords(list): 756 """A class to hold a list of ShapeRecord objects. Subclasses list to ensure compatibility with 757 former work and to reuse all the optimizations of the builtin list. 758 In addition to the list interface, this also provides the GeoJSON __geo_interface__ 759 to return a FeatureCollection dictionary.""" 760 761 def __repr__(self): 762 return 'ShapeRecords: {}'.format(list(self)) 763 764 @property 765 def __geo_interface__(self): 766 return {'type': 'FeatureCollection', 767 'features': [shaperec.__geo_interface__ for shaperec in self]} 768 769class ShapefileException(Exception): 770 """An exception to handle shapefile specific problems.""" 771 pass 772 773class Reader(object): 774 """Reads the three files of a shapefile as a unit or 775 separately. If one of the three files (.shp, .shx, 776 .dbf) is missing no exception is thrown until you try 777 to call a method that depends on that particular file. 778 The .shx index file is used if available for efficiency 779 but is not required to read the geometry from the .shp 780 file. The "shapefile" argument in the constructor is the 781 name of the file you want to open. 782 783 You can instantiate a Reader without specifying a shapefile 784 and then specify one later with the load() method. 785 786 Only the shapefile headers are read upon loading. Content 787 within each file is only accessed when required and as 788 efficiently as possible. Shapefiles are usually not large 789 but they can be. 790 """ 791 def __init__(self, *args, **kwargs): 792 self.shp = None 793 self.shx = None 794 self.dbf = None 795 self.shapeName = "Not specified" 796 self._offsets = [] 797 self.shpLength = None 798 self.numRecords = None 799 self.numShapes = None 800 self.fields = [] 801 self.__dbfHdrLength = 0 802 self.__fieldposition_lookup = {} 803 self.encoding = kwargs.pop('encoding', 'utf-8') 804 self.encodingErrors = kwargs.pop('encodingErrors', 'strict') 805 # See if a shapefile name was passed as the first argument 806 if len(args) > 0: 807 if is_string(args[0]): 808 self.load(args[0]) 809 return 810 # Otherwise, load from separate shp/shx/dbf args (must be file-like) 811 if "shp" in kwargs.keys(): 812 if hasattr(kwargs["shp"], "read"): 813 self.shp = kwargs["shp"] 814 # Copy if required 815 try: 816 self.shp.seek(0) 817 except (NameError, io.UnsupportedOperation): 818 self.shp = io.BytesIO(self.shp.read()) 819 else: 820 raise ShapefileException('The shp arg must be file-like.') 821 822 if "shx" in kwargs.keys(): 823 if hasattr(kwargs["shx"], "read"): 824 self.shx = kwargs["shx"] 825 # Copy if required 826 try: 827 self.shx.seek(0) 828 except (NameError, io.UnsupportedOperation): 829 self.shx = io.BytesIO(self.shx.read()) 830 else: 831 raise ShapefileException('The shx arg must be file-like.') 832 833 if "dbf" in kwargs.keys(): 834 if hasattr(kwargs["dbf"], "read"): 835 self.dbf = kwargs["dbf"] 836 # Copy if required 837 try: 838 self.dbf.seek(0) 839 except (NameError, io.UnsupportedOperation): 840 self.dbf = io.BytesIO(self.dbf.read()) 841 else: 842 raise ShapefileException('The dbf arg must be file-like.') 843 844 # Load the files 845 if self.shp or self.dbf: 846 self.load() 847 848 def __str__(self): 849 """ 850 Use some general info on the shapefile as __str__ 851 """ 852 info = ['shapefile Reader'] 853 if self.shp: 854 info.append(" {} shapes (type '{}')".format( 855 len(self), SHAPETYPE_LOOKUP[self.shapeType])) 856 if self.dbf: 857 info.append(' {} records ({} fields)'.format( 858 len(self), len(self.fields))) 859 return '\n'.join(info) 860 861 def __enter__(self): 862 """ 863 Enter phase of context manager. 864 """ 865 return self 866 867 def __exit__(self, exc_type, exc_val, exc_tb): 868 """ 869 Exit phase of context manager, close opened files. 870 """ 871 self.close() 872 873 def __len__(self): 874 """Returns the number of shapes/records in the shapefile.""" 875 if self.dbf: 876 # Preferably use dbf record count 877 if self.numRecords is None: 878 self.__dbfHeader() 879 880 return self.numRecords 881 882 elif self.shp: 883 # Otherwise use shape count 884 if self.shx: 885 # Use index file to get total count 886 if self.numShapes is None: 887 # File length (16-bit word * 2 = bytes) - header length 888 self.shx.seek(24) 889 shxRecordLength = (unpack(">i", self.shx.read(4))[0] * 2) - 100 890 self.numShapes = shxRecordLength // 8 891 892 return self.numShapes 893 894 else: 895 # Index file not available, iterate all shapes to get total count 896 if self.numShapes is None: 897 i = 0 898 for i,shape in enumerate(self.iterShapes()): 899 i += 1 900 self.numShapes = i 901 902 return self.numShapes 903 904 else: 905 # No file loaded yet, treat as 'empty' shapefile 906 return 0 907 908 def __iter__(self): 909 """Iterates through the shapes/records in the shapefile.""" 910 for shaperec in self.iterShapeRecords(): 911 yield shaperec 912 913 @property 914 def __geo_interface__(self): 915 shaperecords = self.shapeRecords() 916 fcollection = shaperecords.__geo_interface__ 917 fcollection['bbox'] = list(self.bbox) 918 return fcollection 919 920 @property 921 def shapeTypeName(self): 922 return SHAPETYPE_LOOKUP[self.shapeType] 923 924 def load(self, shapefile=None): 925 """Opens a shapefile from a filename or file-like 926 object. Normally this method would be called by the 927 constructor with the file name as an argument.""" 928 if shapefile: 929 (shapeName, ext) = os.path.splitext(shapefile) 930 self.shapeName = shapeName 931 self.load_shp(shapeName) 932 self.load_shx(shapeName) 933 self.load_dbf(shapeName) 934 if not (self.shp or self.dbf): 935 raise ShapefileException("Unable to open %s.dbf or %s.shp." % (shapeName, shapeName)) 936 if self.shp: 937 self.__shpHeader() 938 if self.dbf: 939 self.__dbfHeader() 940 941 def load_shp(self, shapefile_name): 942 """ 943 Attempts to load file with .shp extension as both lower and upper case 944 """ 945 shp_ext = 'shp' 946 try: 947 self.shp = open("%s.%s" % (shapefile_name, shp_ext), "rb") 948 except IOError: 949 try: 950 self.shp = open("%s.%s" % (shapefile_name, shp_ext.upper()), "rb") 951 except IOError: 952 pass 953 954 def load_shx(self, shapefile_name): 955 """ 956 Attempts to load file with .shx extension as both lower and upper case 957 """ 958 shx_ext = 'shx' 959 try: 960 self.shx = open("%s.%s" % (shapefile_name, shx_ext), "rb") 961 except IOError: 962 try: 963 self.shx = open("%s.%s" % (shapefile_name, shx_ext.upper()), "rb") 964 except IOError: 965 pass 966 967 def load_dbf(self, shapefile_name): 968 """ 969 Attempts to load file with .dbf extension as both lower and upper case 970 """ 971 dbf_ext = 'dbf' 972 try: 973 self.dbf = open("%s.%s" % (shapefile_name, dbf_ext), "rb") 974 except IOError: 975 try: 976 self.dbf = open("%s.%s" % (shapefile_name, dbf_ext.upper()), "rb") 977 except IOError: 978 pass 979 980 def __del__(self): 981 self.close() 982 983 def close(self): 984 for attribute in (self.shp, self.shx, self.dbf): 985 if hasattr(attribute, 'close'): 986 try: 987 attribute.close() 988 except IOError: 989 pass 990 991 def __getFileObj(self, f): 992 """Checks to see if the requested shapefile file object is 993 available. If not a ShapefileException is raised.""" 994 if not f: 995 raise ShapefileException("Shapefile Reader requires a shapefile or file-like object.") 996 if self.shp and self.shpLength is None: 997 self.load() 998 if self.dbf and len(self.fields) == 0: 999 self.load() 1000 return f 1001 1002 def __restrictIndex(self, i): 1003 """Provides list-like handling of a record index with a clearer 1004 error message if the index is out of bounds.""" 1005 if self.numRecords: 1006 rmax = self.numRecords - 1 1007 if abs(i) > rmax: 1008 raise IndexError("Shape or Record index out of range.") 1009 if i < 0: i = range(self.numRecords)[i] 1010 return i 1011 1012 def __shpHeader(self): 1013 """Reads the header information from a .shp or .shx file.""" 1014 if not self.shp: 1015 raise ShapefileException("Shapefile Reader requires a shapefile or file-like object. (no shp file found") 1016 shp = self.shp 1017 # File length (16-bit word * 2 = bytes) 1018 shp.seek(24) 1019 self.shpLength = unpack(">i", shp.read(4))[0] * 2 1020 # Shape type 1021 shp.seek(32) 1022 self.shapeType= unpack("<i", shp.read(4))[0] 1023 # The shapefile's bounding box (lower left, upper right) 1024 self.bbox = _Array('d', unpack("<4d", shp.read(32))) 1025 # Elevation 1026 self.zbox = _Array('d', unpack("<2d", shp.read(16))) 1027 # Measure 1028 self.mbox = [] 1029 for m in _Array('d', unpack("<2d", shp.read(16))): 1030 # Measure values less than -10e38 are nodata values according to the spec 1031 if m > NODATA: 1032 self.mbox.append(m) 1033 else: 1034 self.mbox.append(None) 1035 1036 def __shape(self): 1037 """Returns the header info and geometry for a single shape.""" 1038 f = self.__getFileObj(self.shp) 1039 record = Shape() 1040 nParts = nPoints = zmin = zmax = mmin = mmax = None 1041 (recNum, recLength) = unpack(">2i", f.read(8)) 1042 # Determine the start of the next record 1043 next = f.tell() + (2 * recLength) 1044 shapeType = unpack("<i", f.read(4))[0] 1045 record.shapeType = shapeType 1046 # For Null shapes create an empty points list for consistency 1047 if shapeType == 0: 1048 record.points = [] 1049 # All shape types capable of having a bounding box 1050 elif shapeType in (3,5,8,13,15,18,23,25,28,31): 1051 record.bbox = _Array('d', unpack("<4d", f.read(32))) 1052 # Shape types with parts 1053 if shapeType in (3,5,13,15,23,25,31): 1054 nParts = unpack("<i", f.read(4))[0] 1055 # Shape types with points 1056 if shapeType in (3,5,8,13,15,18,23,25,28,31): 1057 nPoints = unpack("<i", f.read(4))[0] 1058 # Read parts 1059 if nParts: 1060 record.parts = _Array('i', unpack("<%si" % nParts, f.read(nParts * 4))) 1061 # Read part types for Multipatch - 31 1062 if shapeType == 31: 1063 record.partTypes = _Array('i', unpack("<%si" % nParts, f.read(nParts * 4))) 1064 # Read points - produces a list of [x,y] values 1065 if nPoints: 1066 flat = unpack("<%sd" % (2 * nPoints), f.read(16*nPoints)) 1067 record.points = list(izip(*(iter(flat),) * 2)) 1068 # Read z extremes and values 1069 if shapeType in (13,15,18,31): 1070 (zmin, zmax) = unpack("<2d", f.read(16)) 1071 record.z = _Array('d', unpack("<%sd" % nPoints, f.read(nPoints * 8))) 1072 # Read m extremes and values 1073 if shapeType in (13,15,18,23,25,28,31): 1074 if next - f.tell() >= 16: 1075 (mmin, mmax) = unpack("<2d", f.read(16)) 1076 # Measure values less than -10e38 are nodata values according to the spec 1077 if next - f.tell() >= nPoints * 8: 1078 record.m = [] 1079 for m in _Array('d', unpack("<%sd" % nPoints, f.read(nPoints * 8))): 1080 if m > NODATA: 1081 record.m.append(m) 1082 else: 1083 record.m.append(None) 1084 else: 1085 record.m = [None for _ in range(nPoints)] 1086 # Read a single point 1087 if shapeType in (1,11,21): 1088 record.points = [_Array('d', unpack("<2d", f.read(16)))] 1089 # Read a single Z value 1090 if shapeType == 11: 1091 record.z = list(unpack("<d", f.read(8))) 1092 # Read a single M value 1093 if shapeType in (21,11): 1094 if next - f.tell() >= 8: 1095 (m,) = unpack("<d", f.read(8)) 1096 else: 1097 m = NODATA 1098 # Measure values less than -10e38 are nodata values according to the spec 1099 if m > NODATA: 1100 record.m = [m] 1101 else: 1102 record.m = [None] 1103 # Seek to the end of this record as defined by the record header because 1104 # the shapefile spec doesn't require the actual content to meet the header 1105 # definition. Probably allowed for lazy feature deletion. 1106 f.seek(next) 1107 return record 1108 1109 def __shapeIndex(self, i=None): 1110 """Returns the offset in a .shp file for a shape based on information 1111 in the .shx index file.""" 1112 shx = self.shx 1113 if not shx: 1114 return None 1115 if not self._offsets: 1116 if self.numShapes is None: 1117 # File length (16-bit word * 2 = bytes) - header length 1118 shx.seek(24) 1119 shxRecordLength = (unpack(">i", shx.read(4))[0] * 2) - 100 1120 self.numShapes = shxRecordLength // 8 1121 # Jump to the first record. 1122 shx.seek(100) 1123 shxRecords = _Array('i') 1124 # Each offset consists of two nrs, only the first one matters 1125 shxRecords.fromfile(shx, 2 * self.numShapes) 1126 if sys.byteorder != 'big': 1127 shxRecords.byteswap() 1128 self._offsets = [2 * el for el in shxRecords[::2]] 1129 if not i == None: 1130 return self._offsets[i] 1131 1132 def shape(self, i=0): 1133 """Returns a shape object for a shape in the geometry 1134 record file.""" 1135 shp = self.__getFileObj(self.shp) 1136 i = self.__restrictIndex(i) 1137 offset = self.__shapeIndex(i) 1138 if not offset: 1139 # Shx index not available so iterate the full list. 1140 for j,k in enumerate(self.iterShapes()): 1141 if j == i: 1142 return k 1143 shp.seek(offset) 1144 return self.__shape() 1145 1146 def shapes(self): 1147 """Returns all shapes in a shapefile.""" 1148 shp = self.__getFileObj(self.shp) 1149 # Found shapefiles which report incorrect 1150 # shp file length in the header. Can't trust 1151 # that so we seek to the end of the file 1152 # and figure it out. 1153 shp.seek(0,2) 1154 self.shpLength = shp.tell() 1155 shp.seek(100) 1156 shapes = Shapes() 1157 while shp.tell() < self.shpLength: 1158 shapes.append(self.__shape()) 1159 return shapes 1160 1161 def iterShapes(self): 1162 """Returns a generator of shapes in a shapefile. Useful 1163 for handling large shapefiles.""" 1164 shp = self.__getFileObj(self.shp) 1165 shp.seek(0,2) 1166 self.shpLength = shp.tell() 1167 shp.seek(100) 1168 while shp.tell() < self.shpLength: 1169 yield self.__shape() 1170 1171 def __dbfHeader(self): 1172 """Reads a dbf header. Xbase-related code borrows heavily from ActiveState Python Cookbook Recipe 362715 by Raymond Hettinger""" 1173 if not self.dbf: 1174 raise ShapefileException("Shapefile Reader requires a shapefile or file-like object. (no dbf file found)") 1175 dbf = self.dbf 1176 # read relevant header parts 1177 self.numRecords, self.__dbfHdrLength, self.__recordLength = \ 1178 unpack("<xxxxLHH20x", dbf.read(32)) 1179 # read fields 1180 numFields = (self.__dbfHdrLength - 33) // 32 1181 for field in range(numFields): 1182 fieldDesc = list(unpack("<11sc4xBB14x", dbf.read(32))) 1183 name = 0 1184 idx = 0 1185 if b"\x00" in fieldDesc[name]: 1186 idx = fieldDesc[name].index(b"\x00") 1187 else: 1188 idx = len(fieldDesc[name]) - 1 1189 fieldDesc[name] = fieldDesc[name][:idx] 1190 fieldDesc[name] = u(fieldDesc[name], self.encoding, self.encodingErrors) 1191 fieldDesc[name] = fieldDesc[name].lstrip() 1192 fieldDesc[1] = u(fieldDesc[1], 'ascii') 1193 self.fields.append(fieldDesc) 1194 terminator = dbf.read(1) 1195 if terminator != b"\r": 1196 raise ShapefileException("Shapefile dbf header lacks expected terminator. (likely corrupt?)") 1197 self.fields.insert(0, ('DeletionFlag', 'C', 1, 0)) 1198 fmt,fmtSize = self.__recordFmt() 1199 self.__recStruct = Struct(fmt) 1200 1201 # Store the field positions 1202 self.__fieldposition_lookup = dict((f[0], i) for i, f in enumerate(self.fields[1:])) 1203 1204 def __recordFmt(self): 1205 """Calculates the format and size of a .dbf record.""" 1206 if self.numRecords is None: 1207 self.__dbfHeader() 1208 fmt = ''.join(['%ds' % fieldinfo[2] for fieldinfo in self.fields]) 1209 fmtSize = calcsize(fmt) 1210 # total size of fields should add up to recordlength from the header 1211 while fmtSize < self.__recordLength: 1212 # if not, pad byte until reaches recordlength 1213 fmt += "x" 1214 fmtSize += 1 1215 return (fmt, fmtSize) 1216 1217 def __record(self, oid=None): 1218 """Reads and returns a dbf record row as a list of values.""" 1219 f = self.__getFileObj(self.dbf) 1220 recordContents = self.__recStruct.unpack(f.read(self.__recStruct.size)) 1221 if recordContents[0] != b' ': 1222 # deleted record 1223 return None 1224 record = [] 1225 for (name, typ, size, deci), value in zip(self.fields, recordContents): 1226 if name == 'DeletionFlag': 1227 continue 1228 elif typ in ("N","F"): 1229 # numeric or float: number stored as a string, right justified, and padded with blanks to the width of the field. 1230 value = value.split(b'\0')[0] 1231 value = value.replace(b'*', b'') # QGIS NULL is all '*' chars 1232 if value == b'': 1233 value = None 1234 elif deci: 1235 try: 1236 value = float(value) 1237 except ValueError: 1238 #not parseable as float, set to None 1239 value = None 1240 else: 1241 # force to int 1242 try: 1243 # first try to force directly to int. 1244 # forcing a large int to float and back to int 1245 # will lose information and result in wrong nr. 1246 value = int(value) 1247 except ValueError: 1248 # forcing directly to int failed, so was probably a float. 1249 try: 1250 value = int(float(value)) 1251 except ValueError: 1252 #not parseable as int, set to None 1253 value = None 1254 elif typ == 'D': 1255 # date: 8 bytes - date stored as a string in the format YYYYMMDD. 1256 if not value.replace(b'\x00', b'').replace(b' ', b'').replace(b'0', b''): 1257 # dbf date field has no official null value 1258 # but can check for all hex null-chars, all spaces, or all 0s (QGIS null) 1259 value = None 1260 else: 1261 try: 1262 # return as python date object 1263 y, m, d = int(value[:4]), int(value[4:6]), int(value[6:8]) 1264 value = date(y, m, d) 1265 except: 1266 # if invalid date, just return as unicode string so user can decide 1267 value = u(value.strip()) 1268 elif typ == 'L': 1269 # logical: 1 byte - initialized to 0x20 (space) otherwise T or F. 1270 if value == b" ": 1271 value = None # space means missing or not yet set 1272 else: 1273 if value in b'YyTt1': 1274 value = True 1275 elif value in b'NnFf0': 1276 value = False 1277 else: 1278 value = None # unknown value is set to missing 1279 else: 1280 # anything else is forced to string/unicode 1281 value = u(value, self.encoding, self.encodingErrors) 1282 value = value.strip() 1283 record.append(value) 1284 1285 return _Record(self.__fieldposition_lookup, record, oid) 1286 1287 def record(self, i=0): 1288 """Returns a specific dbf record based on the supplied index.""" 1289 f = self.__getFileObj(self.dbf) 1290 if self.numRecords is None: 1291 self.__dbfHeader() 1292 i = self.__restrictIndex(i) 1293 recSize = self.__recStruct.size 1294 f.seek(0) 1295 f.seek(self.__dbfHdrLength + (i * recSize)) 1296 return self.__record(oid=i) 1297 1298 def records(self): 1299 """Returns all records in a dbf file.""" 1300 if self.numRecords is None: 1301 self.__dbfHeader() 1302 records = [] 1303 f = self.__getFileObj(self.dbf) 1304 f.seek(self.__dbfHdrLength) 1305 for i in range(self.numRecords): 1306 r = self.__record(oid=i) 1307 if r: 1308 records.append(r) 1309 return records 1310 1311 def iterRecords(self): 1312 """Returns a generator of records in a dbf file. 1313 Useful for large shapefiles or dbf files.""" 1314 if self.numRecords is None: 1315 self.__dbfHeader() 1316 f = self.__getFileObj(self.dbf) 1317 f.seek(self.__dbfHdrLength) 1318 for i in xrange(self.numRecords): 1319 r = self.__record() 1320 if r: 1321 yield r 1322 1323 def shapeRecord(self, i=0): 1324 """Returns a combination geometry and attribute record for the 1325 supplied record index.""" 1326 i = self.__restrictIndex(i) 1327 return ShapeRecord(shape=self.shape(i), record=self.record(i)) 1328 1329 def shapeRecords(self): 1330 """Returns a list of combination geometry/attribute records for 1331 all records in a shapefile.""" 1332 return ShapeRecords(self.iterShapeRecords()) 1333 1334 def iterShapeRecords(self): 1335 """Returns a generator of combination geometry/attribute records for 1336 all records in a shapefile.""" 1337 for shape, record in izip(self.iterShapes(), self.iterRecords()): 1338 yield ShapeRecord(shape=shape, record=record) 1339 1340 1341class Writer(object): 1342 """Provides write support for ESRI Shapefiles.""" 1343 def __init__(self, target=None, shapeType=None, autoBalance=False, **kwargs): 1344 self.target = target 1345 self.autoBalance = autoBalance 1346 self.fields = [] 1347 self.shapeType = shapeType 1348 self.shp = self.shx = self.dbf = None 1349 if target: 1350 if not is_string(target): 1351 raise Exception('The target filepath {} must be of type str/unicode, not {}.'.format(repr(target), type(target)) ) 1352 self.shp = self.__getFileObj(os.path.splitext(target)[0] + '.shp') 1353 self.shx = self.__getFileObj(os.path.splitext(target)[0] + '.shx') 1354 self.dbf = self.__getFileObj(os.path.splitext(target)[0] + '.dbf') 1355 elif kwargs.get('shp') or kwargs.get('shx') or kwargs.get('dbf'): 1356 shp,shx,dbf = kwargs.get('shp'), kwargs.get('shx'), kwargs.get('dbf') 1357 if shp: 1358 self.shp = self.__getFileObj(shp) 1359 if shx: 1360 self.shx = self.__getFileObj(shx) 1361 if dbf: 1362 self.dbf = self.__getFileObj(dbf) 1363 else: 1364 raise Exception('Either the target filepath, or any of shp, shx, or dbf must be set to create a shapefile.') 1365 # Initiate with empty headers, to be finalized upon closing 1366 if self.shp: self.shp.write(b'9'*100) 1367 if self.shx: self.shx.write(b'9'*100) 1368 # Geometry record offsets and lengths for writing shx file. 1369 self.recNum = 0 1370 self.shpNum = 0 1371 self._bbox = None 1372 self._zbox = None 1373 self._mbox = None 1374 # Use deletion flags in dbf? Default is false (0). Note: Currently has no effect, records should NOT contain deletion flags. 1375 self.deletionFlag = 0 1376 # Encoding 1377 self.encoding = kwargs.pop('encoding', 'utf-8') 1378 self.encodingErrors = kwargs.pop('encodingErrors', 'strict') 1379 1380 def __len__(self): 1381 """Returns the current number of features written to the shapefile. 1382 If shapes and records are unbalanced, the length is considered the highest 1383 of the two.""" 1384 return max(self.recNum, self.shpNum) 1385 1386 def __enter__(self): 1387 """ 1388 Enter phase of context manager. 1389 """ 1390 return self 1391 1392 def __exit__(self, exc_type, exc_val, exc_tb): 1393 """ 1394 Exit phase of context manager, finish writing and close the files. 1395 """ 1396 self.close() 1397 1398 def __del__(self): 1399 self.close() 1400 1401 def close(self): 1402 """ 1403 Write final shp, shx, and dbf headers, close opened files. 1404 """ 1405 # Check if any of the files have already been closed 1406 shp_open = self.shp and not (hasattr(self.shp, 'closed') and self.shp.closed) 1407 shx_open = self.shx and not (hasattr(self.shx, 'closed') and self.shx.closed) 1408 dbf_open = self.dbf and not (hasattr(self.dbf, 'closed') and self.dbf.closed) 1409 1410 # Balance if already not balanced 1411 if self.shp and shp_open and self.dbf and dbf_open: 1412 if self.autoBalance: 1413 self.balance() 1414 if self.recNum != self.shpNum: 1415 raise ShapefileException("When saving both the dbf and shp file, " 1416 "the number of records (%s) must correspond " 1417 "with the number of shapes (%s)" % (self.recNum, self.shpNum)) 1418 # Fill in the blank headers 1419 if self.shp and shp_open: 1420 self.__shapefileHeader(self.shp, headerType='shp') 1421 if self.shx and shx_open: 1422 self.__shapefileHeader(self.shx, headerType='shx') 1423 1424 # Update the dbf header with final length etc 1425 if self.dbf and dbf_open: 1426 self.__dbfHeader() 1427 1428 # Close files, if target is a filepath 1429 if self.target: 1430 for attribute in (self.shp, self.shx, self.dbf): 1431 if hasattr(attribute, 'close'): 1432 try: 1433 attribute.close() 1434 except IOError: 1435 pass 1436 1437 def __getFileObj(self, f): 1438 """Safety handler to verify file-like objects""" 1439 if not f: 1440 raise ShapefileException("No file-like object available.") 1441 elif hasattr(f, "write"): 1442 return f 1443 else: 1444 pth = os.path.split(f)[0] 1445 if pth and not os.path.exists(pth): 1446 os.makedirs(pth) 1447 return open(f, "wb+") 1448 1449 def __shpFileLength(self): 1450 """Calculates the file length of the shp file.""" 1451 # Remember starting position 1452 start = self.shp.tell() 1453 # Calculate size of all shapes 1454 self.shp.seek(0,2) 1455 size = self.shp.tell() 1456 # Calculate size as 16-bit words 1457 size //= 2 1458 # Return to start 1459 self.shp.seek(start) 1460 return size 1461 1462 def __bbox(self, s): 1463 x = [] 1464 y = [] 1465 if len(s.points) > 0: 1466 px, py = list(zip(*s.points))[:2] 1467 x.extend(px) 1468 y.extend(py) 1469 else: 1470 # this should not happen. 1471 # any shape that is not null should have at least one point, and only those should be sent here. 1472 # could also mean that earlier code failed to add points to a non-null shape. 1473 raise Exception("Cannot create bbox. Expected a valid shape with at least one point. Got a shape of type '%s' and 0 points." % s.shapeType) 1474 bbox = [min(x), min(y), max(x), max(y)] 1475 # update global 1476 if self._bbox: 1477 # compare with existing 1478 self._bbox = [min(bbox[0],self._bbox[0]), min(bbox[1],self._bbox[1]), max(bbox[2],self._bbox[2]), max(bbox[3],self._bbox[3])] 1479 else: 1480 # first time bbox is being set 1481 self._bbox = bbox 1482 return bbox 1483 1484 def __zbox(self, s): 1485 z = [] 1486 for p in s.points: 1487 try: 1488 z.append(p[2]) 1489 except IndexError: 1490 # point did not have z value 1491 # setting it to 0 is probably ok, since it means all are on the same elavation 1492 z.append(0) 1493 zbox = [min(z), max(z)] 1494 # update global 1495 if self._zbox: 1496 # compare with existing 1497 self._zbox = [min(zbox[0],self._zbox[0]), max(zbox[1],self._zbox[1])] 1498 else: 1499 # first time zbox is being set 1500 self._zbox = zbox 1501 return zbox 1502 1503 def __mbox(self, s): 1504 mpos = 3 if s.shapeType in (11,13,15,18,31) else 2 1505 m = [] 1506 for p in s.points: 1507 try: 1508 if p[mpos] is not None: 1509 # mbox should only be calculated on valid m values 1510 m.append(p[mpos]) 1511 except IndexError: 1512 # point did not have m value so is missing 1513 # mbox should only be calculated on valid m values 1514 pass 1515 if not m: 1516 # only if none of the shapes had m values, should mbox be set to missing m values 1517 m.append(NODATA) 1518 mbox = [min(m), max(m)] 1519 # update global 1520 if self._mbox: 1521 # compare with existing 1522 self._mbox = [min(mbox[0],self._mbox[0]), max(mbox[1],self._mbox[1])] 1523 else: 1524 # first time mbox is being set 1525 self._mbox = mbox 1526 return mbox 1527 1528 @property 1529 def shapeTypeName(self): 1530 return SHAPETYPE_LOOKUP[self.shapeType] 1531 1532 def bbox(self): 1533 """Returns the current bounding box for the shapefile which is 1534 the lower-left and upper-right corners. It does not contain the 1535 elevation or measure extremes.""" 1536 return self._bbox 1537 1538 def zbox(self): 1539 """Returns the current z extremes for the shapefile.""" 1540 return self._zbox 1541 1542 def mbox(self): 1543 """Returns the current m extremes for the shapefile.""" 1544 return self._mbox 1545 1546 def __shapefileHeader(self, fileObj, headerType='shp'): 1547 """Writes the specified header type to the specified file-like object. 1548 Several of the shapefile formats are so similar that a single generic 1549 method to read or write them is warranted.""" 1550 f = self.__getFileObj(fileObj) 1551 f.seek(0) 1552 # File code, Unused bytes 1553 f.write(pack(">6i", 9994,0,0,0,0,0)) 1554 # File length (Bytes / 2 = 16-bit words) 1555 if headerType == 'shp': 1556 f.write(pack(">i", self.__shpFileLength())) 1557 elif headerType == 'shx': 1558 f.write(pack('>i', ((100 + (self.shpNum * 8)) // 2))) 1559 # Version, Shape type 1560 if self.shapeType is None: 1561 self.shapeType = NULL 1562 f.write(pack("<2i", 1000, self.shapeType)) 1563 # The shapefile's bounding box (lower left, upper right) 1564 if self.shapeType != 0: 1565 try: 1566 bbox = self.bbox() 1567 if bbox is None: 1568 # The bbox is initialized with None, so this would mean the shapefile contains no valid geometries. 1569 # In such cases of empty shapefiles, ESRI spec says the bbox values are 'unspecified'. 1570 # Not sure what that means, so for now just setting to 0s, which is the same behavior as in previous versions. 1571 # This would also make sense since the Z and M bounds are similarly set to 0 for non-Z/M type shapefiles. 1572 bbox = [0,0,0,0] 1573 f.write(pack("<4d", *bbox)) 1574 except error: 1575 raise ShapefileException("Failed to write shapefile bounding box. Floats required.") 1576 else: 1577 f.write(pack("<4d", 0,0,0,0)) 1578 # Elevation 1579 if self.shapeType in (11,13,15,18): 1580 # Z values are present in Z type 1581 zbox = self.zbox() 1582 if zbox is None: 1583 # means we have empty shapefile/only null geoms (see commentary on bbox above) 1584 zbox = [0,0] 1585 else: 1586 # As per the ESRI shapefile spec, the zbox for non-Z type shapefiles are set to 0s 1587 zbox = [0,0] 1588 # Measure 1589 if self.shapeType in (11,13,15,18,21,23,25,28,31): 1590 # M values are present in M or Z type 1591 mbox = self.mbox() 1592 if mbox is None: 1593 # means we have empty shapefile/only null geoms (see commentary on bbox above) 1594 mbox = [0,0] 1595 else: 1596 # As per the ESRI shapefile spec, the mbox for non-M type shapefiles are set to 0s 1597 mbox = [0,0] 1598 # Try writing 1599 try: 1600 f.write(pack("<4d", zbox[0], zbox[1], mbox[0], mbox[1])) 1601 except error: 1602 raise ShapefileException("Failed to write shapefile elevation and measure values. Floats required.") 1603 1604 def __dbfHeader(self): 1605 """Writes the dbf header and field descriptors.""" 1606 f = self.__getFileObj(self.dbf) 1607 f.seek(0) 1608 version = 3 1609 year, month, day = time.localtime()[:3] 1610 year -= 1900 1611 # Get all fields, ignoring DeletionFlag if specified 1612 fields = [field for field in self.fields if field[0] != 'DeletionFlag'] 1613 # Ensure has at least one field 1614 if not fields: 1615 raise ShapefileException("Shapefile dbf file must contain at least one field.") 1616 numRecs = self.recNum 1617 numFields = len(fields) 1618 headerLength = numFields * 32 + 33 1619 if headerLength >= 65535: 1620 raise ShapefileException( 1621 "Shapefile dbf header length exceeds maximum length.") 1622 recordLength = sum([int(field[2]) for field in fields]) + 1 1623 header = pack('<BBBBLHH20x', version, year, month, day, numRecs, 1624 headerLength, recordLength) 1625 f.write(header) 1626 # Field descriptors 1627 for field in fields: 1628 name, fieldType, size, decimal = field 1629 name = b(name, self.encoding, self.encodingErrors) 1630 name = name.replace(b' ', b'_') 1631 name = name[:10].ljust(11).replace(b' ', b'\x00') 1632 fieldType = b(fieldType, 'ascii') 1633 size = int(size) 1634 fld = pack('<11sc4xBB14x', name, fieldType, size, decimal) 1635 f.write(fld) 1636 # Terminator 1637 f.write(b'\r') 1638 1639 def shape(self, s): 1640 # Balance if already not balanced 1641 if self.autoBalance and self.recNum < self.shpNum: 1642 self.balance() 1643 # Check is shape or import from geojson 1644 if not isinstance(s, Shape): 1645 if hasattr(s, "__geo_interface__"): 1646 s = s.__geo_interface__ 1647 if isinstance(s, dict): 1648 s = Shape._from_geojson(s) 1649 else: 1650 raise Exception("Can only write Shape objects, GeoJSON dictionaries, " 1651 "or objects with the __geo_interface__, " 1652 "not: %r" % s) 1653 # Write to file 1654 offset,length = self.__shpRecord(s) 1655 self.__shxRecord(offset, length) 1656 1657 def __shpRecord(self, s): 1658 f = self.__getFileObj(self.shp) 1659 offset = f.tell() 1660 # Record number, Content length place holder 1661 self.shpNum += 1 1662 f.write(pack(">2i", self.shpNum, 0)) 1663 start = f.tell() 1664 # Shape Type 1665 if self.shapeType is None and s.shapeType != NULL: 1666 self.shapeType = s.shapeType 1667 if s.shapeType != NULL and s.shapeType != self.shapeType: 1668 raise Exception("The shape's type (%s) must match the type of the shapefile (%s)." % (s.shapeType, self.shapeType)) 1669 f.write(pack("<i", s.shapeType)) 1670 1671 # For point just update bbox of the whole shapefile 1672 if s.shapeType in (1,11,21): 1673 self.__bbox(s) 1674 # All shape types capable of having a bounding box 1675 if s.shapeType in (3,5,8,13,15,18,23,25,28,31): 1676 try: 1677 f.write(pack("<4d", *self.__bbox(s))) 1678 except error: 1679 raise ShapefileException("Failed to write bounding box for record %s. Expected floats." % self.shpNum) 1680 # Shape types with parts 1681 if s.shapeType in (3,5,13,15,23,25,31): 1682 # Number of parts 1683 f.write(pack("<i", len(s.parts))) 1684 # Shape types with multiple points per record 1685 if s.shapeType in (3,5,8,13,15,18,23,25,28,31): 1686 # Number of points 1687 f.write(pack("<i", len(s.points))) 1688 # Write part indexes 1689 if s.shapeType in (3,5,13,15,23,25,31): 1690 for p in s.parts: 1691 f.write(pack("<i", p)) 1692 # Part types for Multipatch (31) 1693 if s.shapeType == 31: 1694 for pt in s.partTypes: 1695 f.write(pack("<i", pt)) 1696 # Write points for multiple-point records 1697 if s.shapeType in (3,5,8,13,15,18,23,25,28,31): 1698 try: 1699 [f.write(pack("<2d", *p[:2])) for p in s.points] 1700 except error: 1701 raise ShapefileException("Failed to write points for record %s. Expected floats." % self.shpNum) 1702 # Write z extremes and values 1703 # Note: missing z values are autoset to 0, but not sure if this is ideal. 1704 if s.shapeType in (13,15,18,31): 1705 try: 1706 f.write(pack("<2d", *self.__zbox(s))) 1707 except error: 1708 raise ShapefileException("Failed to write elevation extremes for record %s. Expected floats." % self.shpNum) 1709 try: 1710 if hasattr(s,"z"): 1711 # if z values are stored in attribute 1712 f.write(pack("<%sd" % len(s.z), *s.z)) 1713 else: 1714 # if z values are stored as 3rd dimension 1715 [f.write(pack("<d", p[2] if len(p) > 2 else 0)) for p in s.points] 1716 except error: 1717 raise ShapefileException("Failed to write elevation values for record %s. Expected floats." % self.shpNum) 1718 # Write m extremes and values 1719 # When reading a file, pyshp converts NODATA m values to None, so here we make sure to convert them back to NODATA 1720 # Note: missing m values are autoset to NODATA. 1721 if s.shapeType in (13,15,18,23,25,28,31): 1722 try: 1723 f.write(pack("<2d", *self.__mbox(s))) 1724 except error: 1725 raise ShapefileException("Failed to write measure extremes for record %s. Expected floats" % self.shpNum) 1726 try: 1727 if hasattr(s,"m"): 1728 # if m values are stored in attribute 1729 f.write(pack("<%sd" % len(s.m), *[m if m is not None else NODATA for m in s.m])) 1730 else: 1731 # if m values are stored as 3rd/4th dimension 1732 # 0-index position of m value is 3 if z type (x,y,z,m), or 2 if m type (x,y,m) 1733 mpos = 3 if s.shapeType in (13,15,18,31) else 2 1734 [f.write(pack("<d", p[mpos] if len(p) > mpos and p[mpos] is not None else NODATA)) for p in s.points] 1735 except error: 1736 raise ShapefileException("Failed to write measure values for record %s. Expected floats" % self.shpNum) 1737 # Write a single point 1738 if s.shapeType in (1,11,21): 1739 try: 1740 f.write(pack("<2d", s.points[0][0], s.points[0][1])) 1741 except error: 1742 raise ShapefileException("Failed to write point for record %s. Expected floats." % self.shpNum) 1743 # Write a single Z value 1744 # Note: missing z values are autoset to 0, but not sure if this is ideal. 1745 if s.shapeType == 11: 1746 # update the global z box 1747 self.__zbox(s) 1748 # then write value 1749 if hasattr(s, "z"): 1750 # if z values are stored in attribute 1751 try: 1752 if not s.z: 1753 s.z = (0,) 1754 f.write(pack("<d", s.z[0])) 1755 except error: 1756 raise ShapefileException("Failed to write elevation value for record %s. Expected floats." % self.shpNum) 1757 else: 1758 # if z values are stored as 3rd dimension 1759 try: 1760 if len(s.points[0]) < 3: 1761 s.points[0].append(0) 1762 f.write(pack("<d", s.points[0][2])) 1763 except error: 1764 raise ShapefileException("Failed to write elevation value for record %s. Expected floats." % self.shpNum) 1765 # Write a single M value 1766 # Note: missing m values are autoset to NODATA. 1767 if s.shapeType in (11,21): 1768 # update the global m box 1769 self.__mbox(s) 1770 # then write value 1771 if hasattr(s, "m"): 1772 # if m values are stored in attribute 1773 try: 1774 if not s.m or s.m[0] is None: 1775 s.m = (NODATA,) 1776 f.write(pack("<1d", s.m[0])) 1777 except error: 1778 raise ShapefileException("Failed to write measure value for record %s. Expected floats." % self.shpNum) 1779 else: 1780 # if m values are stored as 3rd/4th dimension 1781 # 0-index position of m value is 3 if z type (x,y,z,m), or 2 if m type (x,y,m) 1782 try: 1783 mpos = 3 if s.shapeType == 11 else 2 1784 if len(s.points[0]) < mpos+1: 1785 s.points[0].append(NODATA) 1786 elif s.points[0][mpos] is None: 1787 s.points[0][mpos] = NODATA 1788 f.write(pack("<1d", s.points[0][mpos])) 1789 except error: 1790 raise ShapefileException("Failed to write measure value for record %s. Expected floats." % self.shpNum) 1791 # Finalize record length as 16-bit words 1792 finish = f.tell() 1793 length = (finish - start) // 2 1794 # start - 4 bytes is the content length field 1795 f.seek(start-4) 1796 f.write(pack(">i", length)) 1797 f.seek(finish) 1798 return offset,length 1799 1800 def __shxRecord(self, offset, length): 1801 """Writes the shx records.""" 1802 f = self.__getFileObj(self.shx) 1803 try: 1804 f.write(pack(">i", offset // 2)) 1805 except error: 1806 raise ShapefileException('The .shp file has reached its file size limit > 4294967294 bytes (4.29 GB). To fix this, break up your file into multiple smaller ones.') 1807 f.write(pack(">i", length)) 1808 1809 def record(self, *recordList, **recordDict): 1810 """Creates a dbf attribute record. You can submit either a sequence of 1811 field values or keyword arguments of field names and values. Before 1812 adding records you must add fields for the record values using the 1813 field() method. If the record values exceed the number of fields the 1814 extra ones won't be added. In the case of using keyword arguments to specify 1815 field/value pairs only fields matching the already registered fields 1816 will be added.""" 1817 # Balance if already not balanced 1818 if self.autoBalance and self.recNum > self.shpNum: 1819 self.balance() 1820 1821 if recordList: 1822 record = list(recordList) 1823 elif recordDict: 1824 record = [] 1825 for field in self.fields: 1826 if field[0] == 'DeletionFlag': 1827 continue # ignore deletionflag field in case it was specified 1828 if field[0] in recordDict: 1829 val = recordDict[field[0]] 1830 if val is None: 1831 record.append("") 1832 else: 1833 record.append(val) 1834 else: 1835 # Blank fields for empty record 1836 record = ["" for field in self.fields if field[0] != 'DeletionFlag'] 1837 self.__dbfRecord(record) 1838 1839 def __dbfRecord(self, record): 1840 """Writes the dbf records.""" 1841 f = self.__getFileObj(self.dbf) 1842 if self.recNum == 0: 1843 # first records, so all fields should be set 1844 # allowing us to write the dbf header 1845 # cannot change the fields after this point 1846 self.__dbfHeader() 1847 # first byte of the record is deletion flag, always disabled 1848 f.write(b' ') 1849 # begin 1850 self.recNum += 1 1851 fields = (field for field in self.fields if field[0] != 'DeletionFlag') # ignore deletionflag field in case it was specified 1852 for (fieldName, fieldType, size, deci), value in zip(fields, record): 1853 # write 1854 fieldType = fieldType.upper() 1855 size = int(size) 1856 if fieldType in ("N","F"): 1857 # numeric or float: number stored as a string, right justified, and padded with blanks to the width of the field. 1858 if value in MISSING: 1859 value = b"*"*size # QGIS NULL 1860 elif not deci: 1861 # force to int 1862 try: 1863 # first try to force directly to int. 1864 # forcing a large int to float and back to int 1865 # will lose information and result in wrong nr. 1866 value = int(value) 1867 except ValueError: 1868 # forcing directly to int failed, so was probably a float. 1869 value = int(float(value)) 1870 value = format(value, "d")[:size].rjust(size) # caps the size if exceeds the field size 1871 else: 1872 value = float(value) 1873 value = format(value, ".%sf"%deci)[:size].rjust(size) # caps the size if exceeds the field size 1874 elif fieldType == "D": 1875 # date: 8 bytes - date stored as a string in the format YYYYMMDD. 1876 if isinstance(value, date): 1877 value = '{:04d}{:02d}{:02d}'.format(value.year, value.month, value.day) 1878 elif isinstance(value, list) and len(value) == 3: 1879 value = '{:04d}{:02d}{:02d}'.format(*value) 1880 elif value in MISSING: 1881 value = b'0' * 8 # QGIS NULL for date type 1882 elif is_string(value) and len(value) == 8: 1883 pass # value is already a date string 1884 else: 1885 raise ShapefileException("Date values must be either a datetime.date object, a list, a YYYYMMDD string, or a missing value.") 1886 elif fieldType == 'L': 1887 # logical: 1 byte - initialized to 0x20 (space) otherwise T or F. 1888 if value in MISSING: 1889 value = b' ' # missing is set to space 1890 elif value in [True,1]: 1891 value = b'T' 1892 elif value in [False,0]: 1893 value = b'F' 1894 else: 1895 value = b' ' # unknown is set to space 1896 else: 1897 # anything else is forced to string, truncated to the length of the field 1898 value = b(value, self.encoding, self.encodingErrors)[:size].ljust(size) 1899 if not isinstance(value, bytes): 1900 # just in case some of the numeric format() and date strftime() results are still in unicode (Python 3 only) 1901 value = b(value, 'ascii', self.encodingErrors) # should be default ascii encoding 1902 if len(value) != size: 1903 raise ShapefileException( 1904 "Shapefile Writer unable to pack incorrect sized value" 1905 " (size %d) into field '%s' (size %d)." % (len(value), fieldName, size)) 1906 f.write(value) 1907 1908 def balance(self): 1909 """Adds corresponding empty attributes or null geometry records depending 1910 on which type of record was created to make sure all three files 1911 are in synch.""" 1912 while self.recNum > self.shpNum: 1913 self.null() 1914 while self.recNum < self.shpNum: 1915 self.record() 1916 1917 1918 def null(self): 1919 """Creates a null shape.""" 1920 self.shape(Shape(NULL)) 1921 1922 1923 def point(self, x, y): 1924 """Creates a POINT shape.""" 1925 shapeType = POINT 1926 pointShape = Shape(shapeType) 1927 pointShape.points.append([x, y]) 1928 self.shape(pointShape) 1929 1930 def pointm(self, x, y, m=None): 1931 """Creates a POINTM shape. 1932 If the m (measure) value is not set, it defaults to NoData.""" 1933 shapeType = POINTM 1934 pointShape = Shape(shapeType) 1935 pointShape.points.append([x, y, m]) 1936 self.shape(pointShape) 1937 1938 def pointz(self, x, y, z=0, m=None): 1939 """Creates a POINTZ shape. 1940 If the z (elevation) value is not set, it defaults to 0. 1941 If the m (measure) value is not set, it defaults to NoData.""" 1942 shapeType = POINTZ 1943 pointShape = Shape(shapeType) 1944 pointShape.points.append([x, y, z, m]) 1945 self.shape(pointShape) 1946 1947 1948 def multipoint(self, points): 1949 """Creates a MULTIPOINT shape. 1950 Points is a list of xy values.""" 1951 shapeType = MULTIPOINT 1952 points = [points] # nest the points inside a list to be compatible with the generic shapeparts method 1953 self._shapeparts(parts=points, shapeType=shapeType) 1954 1955 def multipointm(self, points): 1956 """Creates a MULTIPOINTM shape. 1957 Points is a list of xym values. 1958 If the m (measure) value is not included, it defaults to None (NoData).""" 1959 shapeType = MULTIPOINTM 1960 points = [points] # nest the points inside a list to be compatible with the generic shapeparts method 1961 self._shapeparts(parts=points, shapeType=shapeType) 1962 1963 def multipointz(self, points): 1964 """Creates a MULTIPOINTZ shape. 1965 Points is a list of xyzm values. 1966 If the z (elevation) value is not included, it defaults to 0. 1967 If the m (measure) value is not included, it defaults to None (NoData).""" 1968 shapeType = MULTIPOINTZ 1969 points = [points] # nest the points inside a list to be compatible with the generic shapeparts method 1970 self._shapeparts(parts=points, shapeType=shapeType) 1971 1972 1973 def line(self, lines): 1974 """Creates a POLYLINE shape. 1975 Lines is a collection of lines, each made up of a list of xy values.""" 1976 shapeType = POLYLINE 1977 self._shapeparts(parts=lines, shapeType=shapeType) 1978 1979 def linem(self, lines): 1980 """Creates a POLYLINEM shape. 1981 Lines is a collection of lines, each made up of a list of xym values. 1982 If the m (measure) value is not included, it defaults to None (NoData).""" 1983 shapeType = POLYLINEM 1984 self._shapeparts(parts=lines, shapeType=shapeType) 1985 1986 def linez(self, lines): 1987 """Creates a POLYLINEZ shape. 1988 Lines is a collection of lines, each made up of a list of xyzm values. 1989 If the z (elevation) value is not included, it defaults to 0. 1990 If the m (measure) value is not included, it defaults to None (NoData).""" 1991 shapeType = POLYLINEZ 1992 self._shapeparts(parts=lines, shapeType=shapeType) 1993 1994 1995 def poly(self, polys): 1996 """Creates a POLYGON shape. 1997 Polys is a collection of polygons, each made up of a list of xy values. 1998 Note that for ordinary polygons the coordinates must run in a clockwise direction. 1999 If some of the polygons are holes, these must run in a counterclockwise direction.""" 2000 shapeType = POLYGON 2001 self._shapeparts(parts=polys, shapeType=shapeType) 2002 2003 def polym(self, polys): 2004 """Creates a POLYGONM shape. 2005 Polys is a collection of polygons, each made up of a list of xym values. 2006 Note that for ordinary polygons the coordinates must run in a clockwise direction. 2007 If some of the polygons are holes, these must run in a counterclockwise direction. 2008 If the m (measure) value is not included, it defaults to None (NoData).""" 2009 shapeType = POLYGONM 2010 self._shapeparts(parts=polys, shapeType=shapeType) 2011 2012 def polyz(self, polys): 2013 """Creates a POLYGONZ shape. 2014 Polys is a collection of polygons, each made up of a list of xyzm values. 2015 Note that for ordinary polygons the coordinates must run in a clockwise direction. 2016 If some of the polygons are holes, these must run in a counterclockwise direction. 2017 If the z (elevation) value is not included, it defaults to 0. 2018 If the m (measure) value is not included, it defaults to None (NoData).""" 2019 shapeType = POLYGONZ 2020 self._shapeparts(parts=polys, shapeType=shapeType) 2021 2022 2023 def multipatch(self, parts, partTypes): 2024 """Creates a MULTIPATCH shape. 2025 Parts is a collection of 3D surface patches, each made up of a list of xyzm values. 2026 PartTypes is a list of types that define each of the surface patches. 2027 The types can be any of the following module constants: TRIANGLE_STRIP, 2028 TRIANGLE_FAN, OUTER_RING, INNER_RING, FIRST_RING, or RING. 2029 If the z (elavation) value is not included, it defaults to 0. 2030 If the m (measure) value is not included, it defaults to None (NoData).""" 2031 shapeType = MULTIPATCH 2032 polyShape = Shape(shapeType) 2033 polyShape.parts = [] 2034 polyShape.points = [] 2035 for part in parts: 2036 # set part index position 2037 polyShape.parts.append(len(polyShape.points)) 2038 # add points 2039 for point in part: 2040 # Ensure point is list 2041 if not isinstance(point, list): 2042 point = list(point) 2043 polyShape.points.append(point) 2044 polyShape.partTypes = partTypes 2045 # write the shape 2046 self.shape(polyShape) 2047 2048 2049 def _shapeparts(self, parts, shapeType): 2050 """Internal method for adding a shape that has multiple collections of points (parts): 2051 lines, polygons, and multipoint shapes. 2052 """ 2053 polyShape = Shape(shapeType) 2054 polyShape.parts = [] 2055 polyShape.points = [] 2056 # Make sure polygon rings (parts) are closed 2057 if shapeType in (5,15,25,31): 2058 for part in parts: 2059 if part[0] != part[-1]: 2060 part.append(part[0]) 2061 # Add points and part indexes 2062 for part in parts: 2063 # set part index position 2064 polyShape.parts.append(len(polyShape.points)) 2065 # add points 2066 for point in part: 2067 # Ensure point is list 2068 if not isinstance(point, list): 2069 point = list(point) 2070 polyShape.points.append(point) 2071 # write the shape 2072 self.shape(polyShape) 2073 2074 def field(self, name, fieldType="C", size="50", decimal=0): 2075 """Adds a dbf field descriptor to the shapefile.""" 2076 if fieldType == "D": 2077 size = "8" 2078 decimal = 0 2079 elif fieldType == "L": 2080 size = "1" 2081 decimal = 0 2082 if len(self.fields) >= 2046: 2083 raise ShapefileException( 2084 "Shapefile Writer reached maximum number of fields: 2046.") 2085 self.fields.append((name, fieldType, size, decimal)) 2086 2087## def saveShp(self, target): 2088## """Save an shp file.""" 2089## if not hasattr(target, "write"): 2090## target = os.path.splitext(target)[0] + '.shp' 2091## self.shp = self.__getFileObj(target) 2092## self.__shapefileHeader(self.shp, headerType='shp') 2093## self.shp.seek(100) 2094## self._shp.seek(0) 2095## chunk = True 2096## while chunk: 2097## chunk = self._shp.read(self.bufsize) 2098## self.shp.write(chunk) 2099## 2100## def saveShx(self, target): 2101## """Save an shx file.""" 2102## if not hasattr(target, "write"): 2103## target = os.path.splitext(target)[0] + '.shx' 2104## self.shx = self.__getFileObj(target) 2105## self.__shapefileHeader(self.shx, headerType='shx') 2106## self.shx.seek(100) 2107## self._shx.seek(0) 2108## chunk = True 2109## while chunk: 2110## chunk = self._shx.read(self.bufsize) 2111## self.shx.write(chunk) 2112## 2113## def saveDbf(self, target): 2114## """Save a dbf file.""" 2115## if not hasattr(target, "write"): 2116## target = os.path.splitext(target)[0] + '.dbf' 2117## self.dbf = self.__getFileObj(target) 2118## self.__dbfHeader() # writes to .dbf 2119## self._dbf.seek(0) 2120## chunk = True 2121## while chunk: 2122## chunk = self._dbf.read(self.bufsize) 2123## self.dbf.write(chunk) 2124 2125## def save(self, target=None, shp=None, shx=None, dbf=None): 2126## """Save the shapefile data to three files or 2127## three file-like objects. SHP and DBF files can also 2128## be written exclusively using saveShp, saveShx, and saveDbf respectively. 2129## If target is specified but not shp, shx, or dbf then the target path and 2130## file name are used. If no options or specified, a unique base file name 2131## is generated to save the files and the base file name is returned as a 2132## string. 2133## """ 2134## # Balance if already not balanced 2135## if shp and dbf: 2136## if self.autoBalance: 2137## self.balance() 2138## if self.recNum != self.shpNum: 2139## raise ShapefileException("When saving both the dbf and shp file, " 2140## "the number of records (%s) must correspond " 2141## "with the number of shapes (%s)" % (self.recNum, self.shpNum)) 2142## # Save 2143## if shp: 2144## self.saveShp(shp) 2145## if shx: 2146## self.saveShx(shx) 2147## if dbf: 2148## self.saveDbf(dbf) 2149## # Create a unique file name if one is not defined 2150## if not shp and not shx and not dbf: 2151## generated = False 2152## if not target: 2153## temp = tempfile.NamedTemporaryFile(prefix="shapefile_",dir=os.getcwd()) 2154## target = temp.name 2155## generated = True 2156## self.saveShp(target) 2157## self.shp.close() 2158## self.saveShx(target) 2159## self.shx.close() 2160## self.saveDbf(target) 2161## self.dbf.close() 2162## if generated: 2163## return target 2164 2165# Begin Testing 2166def test(**kwargs): 2167 import doctest 2168 doctest.NORMALIZE_WHITESPACE = 1 2169 verbosity = kwargs.get('verbose', 0) 2170 if verbosity == 0: 2171 print('Running doctests...') 2172 2173 # ignore py2-3 unicode differences 2174 import re 2175 class Py23DocChecker(doctest.OutputChecker): 2176 def check_output(self, want, got, optionflags): 2177 if sys.version_info[0] == 2: 2178 got = re.sub("u'(.*?)'", "'\\1'", got) 2179 got = re.sub('u"(.*?)"', '"\\1"', got) 2180 res = doctest.OutputChecker.check_output(self, want, got, optionflags) 2181 return res 2182 def summarize(self): 2183 doctest.OutputChecker.summarize(True) 2184 2185 # run tests 2186 runner = doctest.DocTestRunner(checker=Py23DocChecker(), verbose=verbosity) 2187 with open("README.md","rb") as fobj: 2188 test = doctest.DocTestParser().get_doctest(string=fobj.read().decode("utf8").replace('\r\n','\n'), globs={}, name="README", filename="README.md", lineno=0) 2189 failure_count, test_count = runner.run(test) 2190 2191 # print results 2192 if verbosity: 2193 runner.summarize(True) 2194 else: 2195 if failure_count == 0: 2196 print('All test passed successfully') 2197 elif failure_count > 0: 2198 runner.summarize(verbosity) 2199 2200 return failure_count 2201 2202if __name__ == "__main__": 2203 """ 2204 Doctests are contained in the file 'README.md', and are tested using the built-in 2205 testing libraries. 2206 """ 2207 failure_count = test() 2208 sys.exit(failure_count) 2209