iterShapeRecords(fields=...), and use a "
+ "bbox argument so the Reader's spatial index "
+ "skips features whose bounding boxes can't possibly match."
+)
+
+for buf in (src_shp, src_shx, src_dbf):
+ buf.seek(0)
+
+dst_shp, dst_shx, dst_dbf = io.BytesIO(), io.BytesIO(), io.BytesIO()
+
+# Region of interest: a bounding box covering the eastern half.
+region_bbox = [3, -5, 11, 10]
+keep_fields = ["name", "basin", "length_km"]
+
+reader = shapefile.Reader(shp=src_shp, shx=src_shx, dbf=src_dbf)
+writer = shapefile.Writer(shp=dst_shp, shx=dst_shx, dbf=dst_dbf)
+
+# Copy only the fields we want to keep. PyShp returns each field as a
+# plain list, [name, type, size, decimal], not an object with named
+# attributes -- so we index by position. The first entry is always a
+# DeletionFlag, which we skip.
+for field in reader.fields[1:]:
+ if field[0] in keep_fields:
+ writer.field(*field)
+
+kept = []
+for shape_record in reader.iterShapeRecords(
+ bbox=region_bbox, fields=keep_fields,
+):
+ record = shape_record.record
+ if record["basin"] == "North" and record["length_km"] >= 4.0:
+ writer.shape(shape_record.shape)
+ writer.record(*[record[name] for name in keep_fields])
+ kept.append(record["name"])
+
+writer.close()
+reader.close()
+
+note(f"Kept {len(kept)} feature(s): " + ", ".join(kept))
+
+# Read the new shapefile back and visualize the result.
+for buf in (dst_shp, dst_shx, dst_dbf):
+ buf.seek(0)
+
+result = shapefile.Reader(shp=dst_shp, shx=dst_shx, dbf=dst_dbf)
+
+note(
+ f"Output shapefile: {result.shapeTypeName}, "
+ f"{len(result)} feature(s), "
+ f"{len(result.fields) - 1} field(s) "
+ "(notice the notes field was dropped)."
+)
+
+fig, ax = plt.subplots(figsize=(7, 4))
+
+# Draw the original rivers in light gray for context.
+src_shp.seek(0); src_shx.seek(0); src_dbf.seek(0)
+context = shapefile.Reader(shp=src_shp, shx=src_shx, dbf=src_dbf)
+for shape in context.iterShapes():
+ xs, ys = zip(*shape.points)
+ ax.plot(xs, ys, color="lightgray", linewidth=1)
+context.close()
+
+# Overlay the surviving features in bold.
+for shape_record in result.iterShapeRecords():
+ xs, ys = zip(*shape_record.shape.points)
+ ax.plot(xs, ys, linewidth=2.5,
+ label=shape_record.record["name"])
+
+# Show the query bounding box.
+x0, y0, x1, y1 = region_bbox
+ax.plot([x0, x1, x1, x0, x0], [y0, y0, y1, y1, y0],
+ linestyle="--", color="darkorange", label="query bbox")
+
+ax.set_aspect("equal")
+ax.set_title("Long rivers in the North basin, within the query box")
+ax.legend(loc="lower right", fontsize=8)
+fig.tight_layout()
+display(fig, append=True)
+
+result.close()
\ No newline at end of file
diff --git a/examples/pyshp/filter_and_edit/config.toml b/examples/pyshp/filter_and_edit/config.toml
new file mode 100644
index 0000000..4b2116f
--- /dev/null
+++ b/examples/pyshp/filter_and_edit/config.toml
@@ -0,0 +1 @@
+packages = ["pyshp", "matplotlib"]
diff --git a/examples/pyshp/filter_and_edit/setup.py b/examples/pyshp/filter_and_edit/setup.py
new file mode 100644
index 0000000..60d91f9
--- /dev/null
+++ b/examples/pyshp/filter_and_edit/setup.py
@@ -0,0 +1,18 @@
+"""Lightweight setup for the third example."""
+import js
+from pyscript import window, HTML, display as _display
+
+js.alert = window.alert
+
+
+def display(*args, **kwargs):
+ return _display(*args, **kwargs, target=__pyscript_display_target__)
+
+
+def heading(text, level=2):
+ display(HTML(f"{text}
"), append=True) + diff --git a/examples/pyshp/order.json b/examples/pyshp/order.json new file mode 100644 index 0000000..3fc1195 --- /dev/null +++ b/examples/pyshp/order.json @@ -0,0 +1,5 @@ +[ + "write_and_read_points", + "polygons_and_geojson", + "filter_and_edit" +] diff --git a/examples/pyshp/polygons_and_geojson/code.py b/examples/pyshp/polygons_and_geojson/code.py new file mode 100644 index 0000000..b4eaaf9 --- /dev/null +++ b/examples/pyshp/polygons_and_geojson/code.py @@ -0,0 +1,120 @@ +# --------------------------------------------------------------------- +# Polygon shapefiles, holes, and GeoJSON via __geo_interface__. +# --------------------------------------------------------------------- + +import shapefile +import io +import matplotlib.pyplot as plt +from matplotlib.path import Path +from matplotlib.patches import PathPatch + +heading("Tiny parks dataset") +note( + "Polygons in shapefiles must be closed (the last point repeats " + "the first). Crucially, the shapefile format has no flag to " + "distinguish an outer ring from a hole — the only signal is " + "winding order: outer rings clockwise, holes counterclockwise. " + "Get this wrong and tools that consume the file (including " + "PyShp's own GeoJSON conversion) will misread your geometry." +) + +shp_buf, shx_buf, dbf_buf = io.BytesIO(), io.BytesIO(), io.BytesIO() + +# A square park with a square pond cut out, plus a triangular plaza. +# Outer rings are clockwise; holes are counterclockwise. +park_outer = [(0, 0), (10, 0), (10, 10), (0, 10), (0, 0)] # clockwise +pond_hole = [(3, 3), (3, 7), (7, 7), (7, 3), (3, 3)] # counterclockwise +plaza = [(15, 2), (17, 8), (20, 2), (15, 2)] # clockwise + +with shapefile.Writer(shp=shp_buf, shx=shx_buf, dbf=dbf_buf) as writer: + writer.field("name", "C", size=30) + writer.field("area_sqm", "N", decimal=1) + + # A single polygon feature with an outer ring and a hole. + writer.poly([park_outer, pond_hole]) + writer.record("Riverside Park", 84.0) + + # A separate polygon feature with no holes. + writer.poly([plaza]) + writer.record("Civic Plaza", 15.0) + +for buf in (shp_buf, shx_buf, dbf_buf): + buf.seek(0) + +reader = shapefile.Reader(shp=shp_buf, shx=shx_buf, dbf=dbf_buf) + +heading("Inspecting parts and points", level=3) +for shape_record in reader.iterShapeRecords(): + shape = shape_record.shape + name = shape_record.record["name"] + note( + f"{name}: shape type " + f"{shape.shapeTypeName}, "
+ f"{len(shape.points)} points across "
+ f"{len(shape.parts)} ring(s) "
+ f"(part start indices: {list(shape.parts)})."
+ )
+
+heading("GeoJSON for free", level=3)
+note(
+ "Every Shape, Record, and Reader implements "
+ "__geo_interface__, so converting to GeoJSON is "
+ "a one-liner that any geospatial tool can consume."
+)
+geojson = reader.__geo_interface__
+display(HTML(
+ f"type: {geojson['type']}\\n"
+ f"features: {len(geojson['features'])}\\n"
+ f"first geometry type: {geojson['features'][0]['geometry']['type']}"
+ ""
+), append=True)
+
+heading("Plotting polygons with real holes", level=3)
+note(
+ "GeoJSON's nested-ring structure encodes holes semantically, but "
+ "matplotlib needs you to translate that into a compound "
+ "Path: outer ring plus inner rings as sub-paths in a "
+ "single PathPatch. This renders holes as genuinely "
+ "transparent — unlike the common shortcut of overpainting with a "
+ "white fill, which breaks the moment you have a non-white "
+ "background or layered features."
+)
+
+# Plot the polygons using their GeoJSON coordinates.
+fig, ax = plt.subplots(figsize=(7, 4))
+colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
+
+for i, feature in enumerate(geojson["features"]):
+ geom = feature["geometry"]
+ name = feature["properties"]["name"]
+ color = colors[i % len(colors)]
+ # GeoJSON nests rings as [[outer, hole, ...]] for Polygons,
+ # and one level deeper for MultiPolygons.
+ polygons = (
+ geom["coordinates"]
+ if geom["type"] == "MultiPolygon"
+ else [geom["coordinates"]]
+ )
+ for rings in polygons:
+ # Build a compound path: outer ring + holes as sub-paths.
+ # Matplotlib renders the holes as genuinely transparent when
+ # the path alternates winding between outer and inner rings.
+ vertices = []
+ codes = []
+ for ring in rings:
+ vertices.extend(ring)
+ codes.append(Path.MOVETO)
+ codes.extend([Path.LINETO] * (len(ring) - 2))
+ codes.append(Path.CLOSEPOLY)
+ path = Path(vertices, codes)
+ patch = PathPatch(path, facecolor=color, alpha=0.4, label=name)
+ ax.add_patch(patch)
+
+ax.set_aspect("equal")
+ax.autoscale_view()
+ax.set_title("Parks polygons (with a pond-shaped hole)")
+ax.legend(loc="upper right")
+fig.tight_layout()
+display(fig, append=True)
+
+reader.close()
\ No newline at end of file
diff --git a/examples/pyshp/polygons_and_geojson/config.toml b/examples/pyshp/polygons_and_geojson/config.toml
new file mode 100644
index 0000000..4b2116f
--- /dev/null
+++ b/examples/pyshp/polygons_and_geojson/config.toml
@@ -0,0 +1 @@
+packages = ["pyshp", "matplotlib"]
diff --git a/examples/pyshp/polygons_and_geojson/setup.py b/examples/pyshp/polygons_and_geojson/setup.py
new file mode 100644
index 0000000..1e05eaa
--- /dev/null
+++ b/examples/pyshp/polygons_and_geojson/setup.py
@@ -0,0 +1,17 @@
+"""Lightweight setup for the second example."""
+import js
+from pyscript import window, HTML, display as _display
+
+js.alert = window.alert
+
+
+def display(*args, **kwargs):
+ return _display(*args, **kwargs, target=__pyscript_display_target__)
+
+
+def heading(text, level=2):
+ display(HTML(f"{text}
"), append=True) diff --git a/examples/pyshp/write_and_read_points/code.py b/examples/pyshp/write_and_read_points/code.py new file mode 100644 index 0000000..09b5130 --- /dev/null +++ b/examples/pyshp/write_and_read_points/code.py @@ -0,0 +1,82 @@ +""" +A first look at PyShp: write a small Point shapefile entirely in +memory, then read it back to verify the geometry and attributes. + +PyShp is the pure-Python reader/writer for the ESRI Shapefile format. +The PyPI distribution is named `pyshp`, but the import name is +`shapefile`. Docs and source: https://github.com/GeospatialPython/pyshp +""" +from IPython.core.display import display, HTML +# pyshp is imported as the `shapefile` module. +import shapefile +import io + + +heading("A handful of lighthouses") +note( + "A shapefile is really three files working together: " + ".shp (geometry), .shx (index), "
+ "and .dbf (attribute table). PyShp can read and "
+ "write each one to a file-like object, which is exactly what "
+ "we want in the browser."
+)
+
+# Three buffers stand in for the three files on disk.
+shp_buffer = io.BytesIO()
+shx_buffer = io.BytesIO()
+dbf_buffer = io.BytesIO()
+
+# A small dataset: famous lighthouses with their (lon, lat) coordinates.
+lighthouses = [
+ ("Eddystone", -4.2636, 50.1922, 1882),
+ ("Fastnet", -9.6033, 51.3833, 1904),
+ ("Cape Hatteras", -75.5290, 35.2503, 1870),
+ ("Tower of Hercules", -8.4063, 43.3863, 100),
+]
+
+# Build the shapefile by streaming records into the Writer.
+with shapefile.Writer(
+ shp=shp_buffer, shx=shx_buffer, dbf=dbf_buffer,
+) as writer:
+ # Define the attribute schema before adding any records.
+ writer.field("name", "C", size=40) # text
+ writer.field("year", "N", decimal=0) # integer year built
+
+ for name, lon, lat, year in lighthouses:
+ writer.point(lon, lat)
+ writer.record(name, year)
+
+note(f"Wrote {len(lighthouses)} point features to in-memory buffers.")
+
+# Reading is symmetric: hand the buffers to a Reader.
+for buf in (shp_buffer, shx_buffer, dbf_buffer):
+ buf.seek(0)
+
+reader = shapefile.Reader(
+ shp=shp_buffer, shx=shx_buffer, dbf=dbf_buffer,
+)
+
+heading("What did we just write?", level=3)
+note(
+ f"Shape type: {reader.shapeTypeName}. "
+ f"Feature count: {len(reader)}. "
+ f"Bounding box: {tuple(round(c, 3) for c in reader.bbox)}."
+)
+
+heading("Iterating shape/record pairs", level=3)
+rows = []
+for shape_record in reader.iterShapeRecords():
+ (lon, lat) = shape_record.shape.points[0]
+ name = shape_record.record["name"]
+ year = shape_record.record["year"]
+ rows.append(f"| Name | Built | " + "Lon, Lat |
|---|
{text}
"), append=True) +