Pages

Friday, August 20, 2021

Connecting MultiPolygon Edges to Lines

This post is another Spatial Thoughts Academy Weekly Challenge solutions. The challenge is to find the edge of the polygon in a set of building footprints whose center point is closest to a street and then connect it with the nearest point on the street.  Shapely and Fiona were of course up to the challenge.  The basic algorithm is to grab the exterior-most ring of the building footprints, break out the edges, determine the midpoint of each edge, measure the distance to each road until the shortest distance is found, then grab the nearest point on the road, and finally draw a line between the edge midpoint and the nearest road point.  Shapely's nearest_points() method is really versatile here because you can throw any two geometries at it even if they are different types.  Fiona is used to read and write geopackages which I'm really liking lately as more an more people are using them finally. Here's the code:

Thursday, August 12, 2021

Extracting MultiPolygons with Holes using Shapely and Fiona

SpatialThoughts.com recently posted a challenge on LinkedIn to extract only building footprints with holes from a city-wide dataset.  

What made the challenge interesting is the building footprints are MultiPolygons instead of simple Polygons which stumped a few people's attempts. But this challenge is very straight-forward Shapely to select the polygons and Fiona to read and write the data.

The correct feature count if you extract all of the building footprints with holes or polygon inner rings is 45.  The challenge author later demonstrated how to do it in QGIS using a series of expressions.  I think I had the only Python solution. One guy posted a successful answer in R.

The challenge data is a city extracted from OSM: https://lnkd.in/g7dVqzMx

Here is my solution:

Thursday, October 10, 2019

Learning Geospatial Analysis with Python, 3rd Ed.

Third Edition is on the shelves! Geospatial concepts, Geo-python universe, and pound-for-pound still the most pure-python and minimal-dependency examples you’ll find anywhere so somebody somewhere out there will still be able to do the math.  I'm giving away a few eBook versions if you're willing to write a review on Amazon. DM me @SpatialPython with your name and email! https://bit.ly/35gm0xl

Saturday, January 12, 2019

PotreeConverter Mac Binary

Potree is the amazing javascript WebGL library that can effortlessly display multi-million-point lidar
point clouds in a browser using a static web page.  In order to do that, you have to use PotreeConverter which creates an efficient octree of your LAS file.  But the MacOS built-in compilers with XCode are usually several generations behind and don't like the PotreeConverter source. I was finally able to get it to compile so I created a GitHub repository with the binary and the compile steps.  You still need to download the "resources" directory from the original PotreeConverter repository which contains the viewer webpage template. You may also have to deal with issue https://github.com/potree/PotreeConverter/issues/281 but it's an easy fix.  I have no idea how portable the binary is, but hopefully the compile steps will save time for others. 

Thursday, November 16, 2017

Creating MultiPoint Shapefiles with PyShp

Pyshp let's you create any type of shapefile.  Normally a point shapefile has one point per record.  But
Photo credit: VancouverMom.ca
you can also have a group of points tied to a single record in a Multipoint shapefile.  To create a Multipoint shapefile, you just use the "poly" method in the Writer. The poly method isn't just for polygons. It can also create polylines and polypoints (i.e. MultiPoints)!

Here's a simple Multipoint shapefile:


import shapefile

# Create our writer as a multipoint shapefile
w = shapefile.Writer(shapefile.MULTIPOINT)

# We'll create a single dbf field
w.field("NAME","C",40)

# Create a single-part, multi-point shape 
# by declaring the shape
# type after the parts/points list
w.poly([[[3,4],[5,6],[7,8],[9,8]]], shapeType=shapefile.MULTIPOINT)

# Create a record for this feature
w.record("Group1")

# Save the multipoint shapefile
w.save("mpoint")
Repeat the poly and record steps to add additional shapes. Add another nested list of points in the first poly method argument to add more parts to the same record