LiDAR Processing Methodology

The pipeline has two stages:

• Stage A (ArcGIS Pro): raw LiDAR → 1 m DSM and DTM with verified coordinate handling (sections 02–04).
• Stage B (Python): the DSM → slope and aspect → per-building zonal aggregation → energy estimate, written with rasterio, geopandas, and NumPy (sections 05–07).

The same Python module runs both study areas; only the input surface changes.

Source data: USGS 3DEP Lidar Point Cloud, project TX Central B1 2017, tile over the West Austin Neighborhood Group: 24,579,588 points, retrieved as a compressed LAZ file via the USGS 3DEP LidarExplorer.

ArcGIS Pro's Create LAS Dataset rejected the compressed .laz directly (Error 000814), so I used Convert LAS to decompress to .las and build the LAS dataset in one step. I then produced two rasters with LAS Dataset To Raster, differing only in which returns are used and how cells are populated.

A misread of horizontal or vertical units would corrupt every downstream slope and area calculation, so I verified the CRS before deriving anything from the surfaces.

• Projected CRS: NAD83(2011) / UTM Zone 14N (EPSG 6343), linear unit metres.
• Vertical CRS: NAVD88 height, unit metres.
• Implication: both horizontal and vertical are metres, so slope (degrees) and area (m²) compute directly with no unit conversion. EPSG 6343 was carried through to the Python analysis.

The Python stage starts from the 1 m DSM. I derive slope and aspect directly from the elevation grid using the Horn (1981) 3×3 gradient, the same algorithm ArcGIS uses internally. The raster is read with rasterio, nodata is masked to NaN, and the gradient is divided by the true pixel size so slope is in degrees regardless of cell resolution.

# Horn 3x3 gradient: slope & aspect in degrees from any DSM
z = src.read(1).astype("float64"); z[z == src.nodata] = np.nan
dzdx = np.gradient(z, axis=1) / px      # divide by true pixel size (m)
dzdy = np.gradient(z, axis=0) / py
slope  = np.degrees(np.arctan(np.sqrt(dzdx**2 + dzdy**2)))
aspect = np.degrees(np.arctan2(dzdy, -dzdx))   # compass-corrected below

OpenStreetMap building footprints are pulled with osmnx and reprojected to the DSM's CRS with geopandas. For each footprint, rasterio's geometry_mask rasterises the polygon into a boolean mask over the slope and aspect grids, an in-memory zonal-statistics operation. Cells under the slope cap define the usable fraction; the equator-facing aspect score is averaged only over those usable cells.

# For each OSM footprint: mask the raster & reduce to usable stats
mask = geometry_mask([geom], out_shape=slope.shape,
                     transform=transform, invert=True)
s, a = slope[mask], aspect[mask]                 # cells inside this roof
under_cap   = s <= max_slope_deg                  # slope usability filter
frac_usable = under_cap.mean()                  # fraction of roof mountable
aspect_score = _aspect_score(a[under_cap]).mean()  # 1.0 = faces equator

The annual per-building yield uses a first-order model in which every term is a named parameter:

# Energy model: every factor is explicit
usable_area = footprint_area * usable_fraction * frac_usable   # m²
annual_kwh  = (usable_area * annual_ghi_kwh_m2 * aspect_score
               * panel_efficiency * performance_ratio)         # kWh/yr

The Kathmandu transfer adds one Python step the LiDAR side did not need: the open Copernicus DSM arrives in geographic WGS84 (degrees), so it is reprojected to metric UTM 45N with rasterio's calculate_default_transform and reproject before the identical analysis runs. Slope and area are therefore computed in metres for both cities.