Mapping

People were asking how I was extracting the digital terrain model (tree/building elevation) from the NYS GIS LiDAR point clouds (LAS) files

People were asking how I was extracting the digital terrain model (tree/building elevation) from the NYS GIS LiDAR point clouds (LAS) files. While the digital surface models (bare earth) models are widely available in raster digital elevation model GeoTIFF, these files are bare earth and don’t include vegetation or buildings. But you can easily create raster digital elevation model GeoTIFF from the LIDAR point clouds with the PDAL – Point Data Abstraction Library.

Install PDAL

It’s easy if you are running Ubuntu Linux — it’s in standard distribution with 20.04 LTS or later. Most modern distributions of Linux include it in their repositories.

sudo apt install pdal

Download LiDAR Point Clouds

You can download the individual cloud files from NYS GIS FTP site. You can use the LiDAR Shapefile indexes on that website to figure out which file you need. Be aware you may have to download and store a significant amount of data — each roughly 556 acre tile (less then a square mile) is between 500 and 1,000 MB. So you will want to do this somewhere you have a lot of free unlimited data service and storage on your hard drive. And patience while it downloads.

Extract the Digital Terrain Model

You should create a text file called pdal_dtm.json or something similar. The text file should contain:

{ 
	"pipeline": [
		"18TWN220520.las",
		{ 
			"type": "filters.range",
			"limits": "Classification[1:1]"
		},
		{
			"filename":"dsm.tif",
			"gdaldriver":"GTiff",
			"resolution": 0.5,
			"output_type": "max",
			"type":"writers.gdal"
		}
	]
}

This creates the pdal “pipeline” configuration file for the conversion. An explanation of important lines you will need to change:

“18TWN220520.las” – Input point cloud LAS file
“type”: “filters.range” – Tells PDAL to filter points based on request
“limits”: “Classification[1:1]” – Layer to extract *
“filename”:”dsm.tif” – Export file name
“gdaldriver”:”GTiff” – Use GeoTIFF export (same as standard DEM files)
“resolution”: 0.5 – Resolution in meters for export, typically 0.5m or 1m depending on survey
“output_type”: “max” – Highest point reflected back within that point, to get digital surface elevation. You can also use min/max/mean/count/stdev/idw.
“type”:”writers.gdal” – Use the GDAL library for create GeoTIFF

* Layers in NYS GIS Point Clouds – “These point clouds will have at a minimum 2 classifications; Class 1 Unclassified, and Class 2 Bare Earth.” Layer 1 is the digital surface model, it contains building heights and tree/crop cover. Layer 2 “Bare Earth” is the same as what you get from the LiDAR Digital Elevation GeoTIFFs that are widely available on NYSGIS website. Be aware that the Class 2 points may have no-data areas, where no ground elevation was detected due to tree cover or other things blocking ground. The DEM files you download from NYSGIS have these areas filled in with interpolation from surrounding areas.

Once your configuration file is done, run the command:

pdal pipeline pdal_dtm.json

As the Point Cloud files are internally indexed, the export should take only a matter of seconds on a modern desktop computer. The exported DEM/GeoTIFF file will be between 15 to 25 MB, which is much smaller then the point cloud. Then you can load in Quantum GIS or your favorite non-free software GIS client like any other Raster DEM file.

To calculate the building or tree height, just subtract the Class 1 raster from Class 2 raster, using the Raster Calculator in QGIS or your favorite GIS program.Β  Or better yet, subtract Class 1 points from the already processed DEM files on NYSGIS website, as they have the missing bare earth data filled in.

How to deal with Shapefiles missing Projection information in QGIS

How to deal with Shapefiles missing Projection information in QGIS πŸ—Ί

Your best bet is to find the proper projection system by including the prj file in the same directory as the shapefile. If not available, try looking for the projection information (EPSG code) from where you got the file or the included XML meta data which can be opened in your web browser.

As a last resort, you can use QGIS Find Projection tool (use the Type to locate box in the lower left) after loading the shapefile to get potential projections. That won’t get you one answer but will narrow down the options. Be aware that using the wrong datum with the right coordinate system will introduce small errors in the projection – often a couple hundred meters or less, so if things are slightly off try a different datum. But try to find the official projection from the prj or xml file.

Store your own data as geopackages to avoid lost prj files and future projection issues.

How to label a point in QGIS based on a vector layer below it

How to label a point in QGIS based on a vector layer below it

For example, you might have a point layer and want to label what town it is in. You might have the Census Tiger County Subdivision layer named ‘tl_2020_36_cousub20’, with the field ‘NAMELSAD20’ that you want to return from your point or polygon layer.

attribute(array_first(overlay_nearest(‘tl_2020_36_cousub20′,$currentfeature)),’NAMELSAD20’)

This will give you the nearest feature.

You need to wrap it in array_first as shown above, as overlay_nearest returns an array with one element. Then attribute just returns the attribute based on the object provided.

Alternatives include overlay_contains, overlay_crosses, overlay_intersects, overlay_touches, overlay_within which will return what you want. Try searching the QGIS Help in the QGIS Expressions builder.

You can also use this style of expression to create a spatial join or nearest neighbor join, inside the field calculator.

I’ve noticed that QGIS doesn’t always do a good job simplifying geometries when exporting to an SVG file

I’ve noticed that QGIS doesn’t always do a good job simplifying geometries when exporting to an SVG file. While you can manually simplify layers, an easier method is as follows:

1) Open Layer Styling

2) Go to Symbology

3) Switch fill to Geometry Generator

4) Use the Simplify Expression to Simplify on the Fly, such as simplifying by 100 map units, which when projected in UTM would be meters.

simplify($geometry,100)

How I got interested in GIS πŸ—Ί

How I got interested in GIS πŸ—Ί

I have no formal GIS training, as things were still pretty primitive back when I was in college in the early 2000s especially when it came to web services, online data and open source software. Computers where a lot less powerful back then. I remember vaguely hearing a bit about Remote Sensing when I was involved in the Environmental Science Club in college, but it wasn’t something I ever used. 

I picked up QGIS in 2010, as I was looking for a way to make my own topographic maps, as I wasn’t happy with what you could get on the Internet. I learned I could FOIL the primative campsite shapefile from the DEC and get data from there to help find campsites. I was pretty good at map and compass stuff from  my years in Boy Scouts. Over time, I’ve branched out into other GIS areas. I’m always been very interested in land use but also demographics.

More recently I’ve been doing more automation of processes, using Python and R statistical language to do some map plotting and a lot of Census data gathering and processing. I like working with R it’s fast and easy to implement code in. I’ve also lately been doing a lot more with Leaflet and web services more generally. 

I don’t do anything professionally with maps, it’s just a hobby.