Mapping

Creating Digital Surface Models GeoTIFF Using National Map Downloader, LiDAR Point Clouds and PDAL 🗺

|

Creating Digital Surface Models GeoTIFF Using National Map Downloader, LiDAR Point Clouds and PDAL 🗺

Find LiDAR Point Clouds on the National Map Downloader: https://apps.nationalmap.gov/downloader/#/

Select Find Elevation Source Data (3DEP) – Lidar, IfSAR

Search for LiDAR Point Cloud (LPC)

Click Export all results as a TXT file and save to a directory.

Then run this Unix command on the text file to download point clouds:

cat data.txt | tr -d '\n' | xargs -d '\r' -I {} wget -c -nc '{}'

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):

{ 
    "pipeline": [
        { 
            "type": "readers.las"
        },
        {
            "type": "filters.range",
            "limits": "Classification[1:1]"
        },
        {
            "gdaldriver":"GTiff",
            "resolution": 1.0,
            "output_type": "max",
            "type":"writers.gdal"
        }
    ]
}

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:

ls *.laz | xargs -I {} basename {} .laz | xargs -P3 -I {} pdal pipeline pipeline.txt --readers.las.filename={}.laz --writers.gdal.filename={}.tif

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.

gdalbuildvrt dsm.vrt *.tif

Python code, takes a list of addresses, runs them through NYS Address Management system, then uses GeoPandas to calculate the municipality and Assembly district from GeoPackages then writes out to a CSV file

|

Python, the State Geocoder and Outputting Assembly Districts

Python code, takes a list of addresses, runs them through NYS Address Management system, then uses GeoPandas to calculate the municipality and Assembly district from GeoPackages then writes out to a CSV file. Nothing fancy, but it does a job. I could use it with any other shapefile I want. GeoPandas is nice because it’s fast and doesn’t require loading QGIS.

#!/usr/bin/python

import requests,sys,json,os,csv

import pandas as pd
import geopandas as gpd

lines=[]

# read list of addresses
with open(sys.argv[-1], newline='') as csvfile:	
	for line in csv.DictReader(csvfile):
		lines.append(line)
		
# build address query
query = '{"records": ['
i=0
for line in lines:
	query += '{ "attributes": { "OBJECTID":'+str(i)+', "SINGLELINE": "'+line['Address'].rstrip()+'"} },'+"\n"
	i+=1	
query += ']}'

post = { 'f':'pjson', 'outSR': 4326, 'addresses': query }
url = 'https://gisservices.its.ny.gov/arcgis/rest/services/Locators/Street_and_Address_Composite/GeocodeServer/geocodeAddresses'

# send request to state geocoder
req = requests.post(url, data = post)
locations = json.loads(req.text)['locations']

# parse response
for loc in locations:
	i = loc['attributes']['ResultID']
	lines[i]['y'] = loc['location']['y']
	lines[i]['x'] = loc['location']['x']
	lines[i]['Match_addr'] = loc['attributes']['Match_addr']
	
	# hackish, might cause problems but keeps joins from erroring
	if (lines[i]['x'] == 'NaN'):
		lines[i]['x'] = 0
		lines[i]['y'] = 0

# convert to pandas
locPd = pd.DataFrame(lines,columns=lines[0].keys())
locPd = gpd.GeoDataFrame(locPd,  geometry=gpd.points_from_xy(locPd.x.astype('float32'), locPd.y.astype('float32')))

# add county municipality column
cosub = gpd.read_file(r'/home/andy/Documents/GIS.Data/geocode/cosub.gpkg')
locPd = gpd.sjoin(locPd, cosub, op="within")

del locPd['index_right']

# add ads column
ad = gpd.read_file(r'/home/andy/Documents/GIS.Data/geocode/ad.gpkg')
locPd = gpd.sjoin(locPd, ad, op="within")

del locPd['index_right']

# add sd column
sd = gpd.read_file(r'/home/andy/Documents/GIS.Data/geocode/sd.gpkg')
locPd = gpd.sjoin(locPd, sd, op="within")

del locPd['index_right']

# add cd column
sd = gpd.read_file(r'/home/andy/Documents/GIS.Data/geocode/cd.gpkg')
locPd = gpd.sjoin(locPd, sd, 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[-1])[0]+'-output.csv', index = False, header=True)