Another Example of GeoPandas and Python Spatial Joins 🗺
Yesterday, I posted a much more complicated piece of code that pulled addresses from the SAM (State Address Management) database and did a spatial join to add a column to the file with Assembly District and Municipality. This was a bit too complex, so I made a simpler one for other purposes that doesn’t require the coordinates to be obtained from SAM.
This python script takes two parameters:
Path to a CSV file that contains an X and Y coordinate
Path to a Shapefile or Geopackage to Join Against
Then the code will create a new CSV file with the spatially joined attributes pulled from the Shapefile. I have only run it on a few large data sets, but I found it took roughly 1 second to join 1,000 records from call to end of end of script.
#!/usr/bin/python
import requests,sys,json,os,csv
import pandas as pd
import geopandas as gpd
lines=[]
# read list of addresses from parameter 1
with open(sys.argv[-2], newline='') as csvfile:
for line in csv.DictReader(csvfile):
lines.append(line)
# convert to pandas
locPd = pd.DataFrame(lines,columns=lines[0].keys())
locPd.convert_dtypes()
locPd = gpd.GeoDataFrame(locPd, geometry=gpd.points_from_xy(locPd.x.astype('float32'), locPd.y.astype('float32')))
# run spatial joins against parameter 2
ad = gpd.read_file(sys.argv[-1])
locPd = gpd.sjoin(locPd, ad, op="within")
# remove added geometery and index columns
del locPd['geometry']
del locPd['index_right']
# write pandas back to out csv
locPd.to_csv (os.path.splitext(sys.argv[-2])[0]+'-output.csv', index = False, header=True)
Back when I was in Boy Scouts, we often used 1:24k quadrangle topographic maps from the USGS. These maps from the pre computer age, manually drafted, they are no longer updated in favor of GIS data sets and the computer created National Map product. Some of the trails and landmarks not in the GNIS dataset never made it over to the new maps - especially on state land. And the National Map topographic maps are hardly the works of art that many the old topographic maps were.
Next create a pipeline.txt for pdal with the classification (For DSM, 1 are unclassified points like buildings and treetops, while 2 are ground points, if you want a DEM, you can also make them this way too with Classification of 2:2):
Next convert the point clouds into digital surface model (GeoTIFF), you can use this shell command with xargs to go over each LAS file, using the above pipeline:
The above command can be somewhat slow depending on how many LAZ point clouds you need to go through and your selected resolution. -P3 sets the number of parallel process (3) which can help speed things up a bit.
Now we have the digital surface models raster that can be used in QGIS for hillshade for 3D.
Build a virtual raster (dsm.vrt) for easy loading into QGIS rather than loading separate files.
Lately I have been experimenting with digital terrain models extracted from the LIDAR point cloud. It allows me to measure the height of buildings and tree cover and can be used for better 3d modeling.
One should never confuse the Town of Highlands in Orange County (containing West Point, Bear Mountain, and Bear Mountain Bridge on the Hudson River) with the Town of Highland in Sullivan County (containing Barryville, Eldred, Highland Lake, Minisink Ford and Yulan on Delaware River), with the Hamlet of Highland in the Town of Lloyd in Ulster County (which also borders the Hudson River, and contains the Mid-Hudson Bridge to Poughkeespie).
The hamlet of Highland is pronounced as two seperate words — “high land”, stressing the syllables. In contrast, the Town of Highland in Sullivan County and Highlands in Orange County rhymes with “island”.
You can also drive from Ulster County to Orange County to Sullivan County, without ever crossing any other counties.
