|
Login | Get Token |
|
Los Angeles County IMAGE Directory > services > LARIAC7 > TREE_CANOPY_2024 (ImageServer) |
Help | API Reference
EGISAGSIMAGE1 |
| JSON | SOAP |
This layer is a combination of two datasets. One of which is an Normalized Difference Vegetation Index (NDVI) product derived from the LARIAC7 program and the other is a Normalized Digital Surface Model from our County partners at NV5. These two datasets were combined together to produce a new countywide tree canopy layer. The imagery was collected in Summer of 2024 and only includes urban/suburban areas of the County of Los Angeles.
1. Los Angeles Region Imagery Acquisition Consortium (LARIAC-7) Normalized Difference Vegetation Index (NDVI) tiff tiles. NDVI is a simple numerical indicator that can be used to analyze remote sensing measurements, typically but not necessarily from a space platform, and assess whether the target being observed contains live green vegetation or not. NDVI was developed from the 2024 Color Infrared (CIR) imagery at 3-inch resolution which contains the Red and Near Infrared (NIR) Bands. The 3-inch data was resampled to 1 meter resolution. Click here for more information.
2. Normalized Digital Surface Model (NDSM) tiff tiles. NDSM (normalized digital surface model) represents the relative height of features above the surrounding ground surface. It is a valuable elevation product that provides insights into the relative height of features compared to the surrounding ground. In contrast, DSMs (digital surface models) and DTMs (digital terrain models) represent absolute height referenced to a datum, typically the mean sea level.
Processing Steps:
1. To create NDVI, the following code was used outside of ArcPro to efficiently combine the Red and Infrared bands of the CIR imagery.
import os
import time
import numpy as np
import rasterio
from tqdm import tqdm # progress bar
# =========================
# CONFIG SECTION (EDIT ME)
# =========================
RED_RASTER = r"D:\NDVI\SOURCE\NDVI_CALC_2024\NDVI_CALC_2024\Red.tif"
NIR_RASTER = r"D:\NDVI\SOURCE\NDVI_CALC_2024\NDVI_CALC_2024\Infrared.tif"
OUTPUT_NDVI = r"D:\NDVI\SOURCE\NDVI_CALC_2024\NDVI_CALC_2024\NDVI_2024_V2.tif"
NDVI_NODATA = -9999.0
OUTPUT_COMPRESS = "NONE" # FASTEST. Switch to DEFLATE later if needed.
BLOCK_SIZE = 512 # Good block alignment for most rasters.
os.environ.setdefault("GDAL_NUM_THREADS", "ALL_CPUS")
def compute_ndvi_streaming():
t0 = time.time()
print(f"\nOpening RED raster: {RED_RASTER}")
with rasterio.open(RED_RASTER) as red_ds:
w = red_ds.width
h = red_ds.height
red_nodata = red_ds.nodata
print(f" Size: {w} × {h}")
with rasterio.open(NIR_RASTER) as nir_ds:
if (nir_ds.width, nir_ds.height) != (w, h):
raise ValueError("Input rasters differ in size.")
if nir_ds.transform != red_ds.transform:
raise ValueError("Input rasters have different transforms.")
if nir_ds.crs != red_ds.crs:
raise ValueError("Input rasters differ in CRS.")
nir_nodata = nir_ds.nodata
profile = red_ds.profile.copy()
profile.update(
driver="GTiff",
dtype="float32",
count=1,
nodata=NDVI_NODATA,
tiled=True,
blockxsize=BLOCK_SIZE,
blockysize=BLOCK_SIZE,
compress=OUTPUT_COMPRESS,
BIGTIFF="IF_SAFER",
)
print("\nCreating output NDVI raster...")
with rasterio.open(OUTPUT_NDVI, "w", **profile) as out_ds:
# Count blocks for progress bar
block_list = list(red_ds.block_windows(1))
print(f"Total blocks to process: {len(block_list)}")
# Progress bar loop
for ji, window in tqdm(block_list, desc="Computing NDVI", unit="block"):
# Read both windows
red = red_ds.read(1, window=window).astype("float32")
nir = nir_ds.read(1, window=window).astype("float32")
# Nodata mask
mask = np.zeros(red.shape, dtype=bool)
if red_nodata is not None:
mask |= (red == red_nodata)
if nir_nodata is not None:
mask |= (nir == nir_nodata)
denom = nir + red
safe = denom != 0
valid = safe & (~mask)
ndvi_block = np.full(red.shape, NDVI_NODATA, dtype="float32")
ndvi_block[valid] = (nir[valid] - red[valid]) / denom[valid]
# Clamp
np.clip(ndvi_block, -1.0, 1.0, out=ndvi_block)
out_ds.write(ndvi_block, 1, window=window)
total = (time.time() - t0) / 60
print(f"\nFinished! Total time: {total:.2f} minutes")
print(f"NDVI saved to: {OUTPUT_NDVI}")
if __name__ == "__main__":
compute_ndvi_streaming()
2. All areas where the surface height in the NDSM Dataset was greater than 3.048 meters or 10 feet*, and the NDVI value was <.1 were extracted and assigned the value 1. This is the trees layer. The formula is listed below.
SetNull("Mosaic_Height_Float2", SetNull("NDVI_MERGED_Band_1",1,"Value <.1"),"Value < 3.048")
3. To provide improved symbolization of trees across the County, and to remove a number of artifacts, a set of additional processes were run. The MajorityFilter command was run twice on the trees layer. For the options, I used 8 neighbors with the majority replacement threshold.
*Important Note: There is not one single standard threshold used across the board for defining tree height. A height of a 3 meter threshold which is used in municipalities in other countries and in research, has been the standard for tree height here in LA County PW.
See a few resources of research and municipalities using the 10ft / 3meter threshold:
https://www.sciencedirect.com/science/article/pii/S2667393225000018 (section 2.3.1, includes a few references)
https://www.mdpi.com/2072-4292/13/4/767 (section 2.3.2)
Here is more info related to the 2016 tree canopy data for the county. They don't have their methodology publicly available, but you might be able to get more info if you reached out to the email listed at the bottom of the page on the second link. I'm pretty sure they also used a 10ft threshold, but without methodology I can't say for certain.
https://treepeople.org/project/la-treecanopy-data/
https://storymaps.arcgis.com/stories/df083f2adb6a4650a738dbf2805674e2
This layer is a combination of two datasets. One of which is an Normalized Difference Vegetation Index (NDVI) product derived from the LARIAC7 program and the other is a Normalized Digital Surface Model from our County partners at NV5. These two datasets were combined together to produce a new countywide tree canopy layer. The imagery was collected in Summer of 2024 and only includes urban/suburban areas of the County of Los Angeles.
1. Los Angeles Region Imagery Acquisition Consortium (LARIAC-7) Normalized Difference Vegetation Index (NDVI) tiff tiles. NDVI is a simple numerical indicator that can be used to analyze remote sensing measurements, typically but not necessarily from a space platform, and assess whether the target being observed contains live green vegetation or not. NDVI was developed from the 2024 Color Infrared (CIR) imagery at 3-inch resolution which contains the Red and Near Infrared (NIR) Bands. The 3-inch data was resampled to 1 meter resolution. Click here for more information.
2. Normalized Digital Surface Model (NDSM) tiff tiles. NDSM (normalized digital surface model) represents the relative height of features above the surrounding ground surface. It is a valuable elevation product that provides insights into the relative height of features compared to the surrounding ground. In contrast, DSMs (digital surface models) and DTMs (digital terrain models) represent absolute height referenced to a datum, typically the mean sea level.
Processing Steps:
1. To create NDVI, the following code was used outside of ArcPro to efficiently combine the Red and Infrared bands of the CIR imagery.
import os
import time
import numpy as np
import rasterio
from tqdm import tqdm # progress bar
# =========================
# CONFIG SECTION (EDIT ME)
# =========================
RED_RASTER = r"Red.tif"
NIR_RASTER = r"Infrared.tif"
OUTPUT_NDVI = r"NDVI_2024_V2.tif"
NDVI_NODATA = -9999.0
OUTPUT_COMPRESS = "NONE" # FASTEST. Switch to DEFLATE later if needed.
BLOCK_SIZE = 512 # Good block alignment for most rasters.
os.environ.setdefault("GDAL_NUM_THREADS", "ALL_CPUS")
def compute_ndvi_streaming():
t0 = time.time()
print(f"\nOpening RED raster: {RED_RASTER}")
with rasterio.open(RED_RASTER) as red_ds:
w = red_ds.width
h = red_ds.height
red_nodata = red_ds.nodata
print(f" Size: {w} × {h}")
with rasterio.open(NIR_RASTER) as nir_ds:
if (nir_ds.width, nir_ds.height) != (w, h):
raise ValueError("Input rasters differ in size.")
if nir_ds.transform != red_ds.transform:
raise ValueError("Input rasters have different transforms.")
if nir_ds.crs != red_ds.crs:
raise ValueError("Input rasters differ in CRS.")
nir_nodata = nir_ds.nodata
profile = red_ds.profile.copy()
profile.update(
driver="GTiff",
dtype="float32",
count=1,
nodata=NDVI_NODATA,
tiled=True,
blockxsize=BLOCK_SIZE,
blockysize=BLOCK_SIZE,
compress=OUTPUT_COMPRESS,
BIGTIFF="IF_SAFER",
)
print("\nCreating output NDVI raster...")
with rasterio.open(OUTPUT_NDVI, "w", **profile) as out_ds:
# Count blocks for progress bar
block_list = list(red_ds.block_windows(1))
print(f"Total blocks to process: {len(block_list)}")
# Progress bar loop
for ji, window in tqdm(block_list, desc="Computing NDVI", unit="block"):
# Read both windows
red = red_ds.read(1, window=window).astype("float32")
nir = nir_ds.read(1, window=window).astype("float32")
# Nodata mask
mask = np.zeros(red.shape, dtype=bool)
if red_nodata is not None:
mask |= (red == red_nodata)
if nir_nodata is not None:
mask |= (nir == nir_nodata)
denom = nir + red
safe = denom != 0
valid = safe & (~mask)
ndvi_block = np.full(red.shape, NDVI_NODATA, dtype="float32")
ndvi_block[valid] = (nir[valid] - red[valid]) / denom[valid]
# Clamp
np.clip(ndvi_block, -1.0, 1.0, out=ndvi_block)
out_ds.write(ndvi_block, 1, window=window)
total = (time.time() - t0) / 60
print(f"\nFinished! Total time: {total:.2f} minutes")
print(f"NDVI saved to: {OUTPUT_NDVI}")
if __name__ == "__main__":
compute_ndvi_streaming()
2. All areas where the surface height in the NDSM Dataset was greater than 3.048 meters or 10 feet*, and the NDVI value was <.1 were extracted and assigned the value 1. This is the trees layer. The formula is listed below.
SetNull("Mosaic_Height_Float2", SetNull("NDVI_MERGED_Band_1",1,"Value <.1"),"Value < 3.048")
3. To provide improved symbolization of trees across the County, and to remove a number of artifacts, a set of additional processes were run. The MajorityFilter command was run twice on the trees layer. For the options, I used 8 neighbors with the majority replacement threshold.
*Important Note: There is not one single standard threshold used across the board for defining tree height. A height of a 3 meter threshold which is used in municipalities in other countries and in research, has been the standard for tree height here in LA County PW.
See a few resources of research and municipalities using the 10ft / 3meter threshold:
https://www.sciencedirect.com/science/article/pii/S2667393225000018 (section 2.3.1, includes a few references)
https://www.mdpi.com/2072-4292/13/4/767 (section 2.3.2)
Here is more info related to the 2016 tree canopy data for the county. They don't have their methodology publicly available, but you might be able to get more info if you reached out to the email listed at the bottom of the page on the second link. I'm pretty sure they also used a 10ft threshold, but without methodology I can't say for certain.
https://treepeople.org/project/la-treecanopy-data/
https://storymaps.arcgis.com/stories/df083f2adb6a4650a738dbf2805674e2