{
 "currentVersion": 10.91,
 "serviceDescription": "<DIV STYLE=\"text-align:Left;\"><DIV><DIV><P><SPAN>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. <\/SPAN><\/P><P><SPAN /><\/P><P STYLE=\"text-indent:20;\"><SPAN>1. <\/SPAN><A href=\"https://lariac-lacounty.hub.arcgis.com/\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>Los Angeles Region Imagery Acquisition Consortium (LARIAC-7)  <\/SPAN><\/A><SPAN>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. <\/SPAN><A href=\"https://gisgeography.com/ndvi-normalized-difference-vegetation-index/\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>Click here for more information.<\/SPAN><\/A><\/P><P><SPAN /><\/P><P STYLE=\"text-indent:20;\"><SPAN>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.<\/SPAN><\/P><P><SPAN /><\/P><P><SPAN STYLE=\"font-weight:bold;\">Processing Steps:<\/SPAN><\/P><P><SPAN>1. To create NDVI, the following code was used outside of ArcPro to efficiently combine the Red and Infrared bands of the CIR imagery.<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>import os<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>import time<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>import numpy as np<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>import rasterio<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>from tqdm import tqdm   # progress bar<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN># =========================<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN># CONFIG SECTION (EDIT ME)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN># =========================<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>RED_RASTER = r\"D:\\NDVI\\SOURCE\\NDVI_CALC_2024\\NDVI_CALC_2024\\Red.tif\"<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>NIR_RASTER = r\"D:\\NDVI\\SOURCE\\NDVI_CALC_2024\\NDVI_CALC_2024\\Infrared.tif\"<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>OUTPUT_NDVI = r\"D:\\NDVI\\SOURCE\\NDVI_CALC_2024\\NDVI_CALC_2024\\NDVI_2024_V2.tif\"<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>NDVI_NODATA = -9999.0<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>OUTPUT_COMPRESS = \"NONE\"        # FASTEST. Switch to DEFLATE later if needed.<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>BLOCK_SIZE = 512                # Good block alignment for most rasters.<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>os.environ.setdefault(\"GDAL_NUM_THREADS\", \"ALL_CPUS\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>def compute_ndvi_streaming():<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>    t0 = time.time()<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>    print(f\"\\nOpening RED raster: {RED_RASTER}\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>    with rasterio.open(RED_RASTER) as red_ds:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>        w = red_ds.width<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>        h = red_ds.height<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>        red_nodata = red_ds.nodata<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>        print(f\"  Size: {w} × {h}\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>        <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>        with rasterio.open(NIR_RASTER) as nir_ds:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            if (nir_ds.width, nir_ds.height) != (w, h):<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                raise ValueError(\"Input rasters differ in size.\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            if nir_ds.transform != red_ds.transform:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                raise ValueError(\"Input rasters have different transforms.\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            if nir_ds.crs != red_ds.crs:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                raise ValueError(\"Input rasters differ in CRS.\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            nir_nodata = nir_ds.nodata<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            profile = red_ds.profile.copy()<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            profile.update(<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                driver=\"GTiff\",<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                dtype=\"float32\",<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                count=1,<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                nodata=NDVI_NODATA,<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                tiled=True,<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                blockxsize=BLOCK_SIZE,<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                blockysize=BLOCK_SIZE,<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                compress=OUTPUT_COMPRESS,<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                BIGTIFF=\"IF_SAFER\",<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            )<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            print(\"\\nCreating output NDVI raster...\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            with rasterio.open(OUTPUT_NDVI, \"w\", **profile) as out_ds:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                # Count blocks for progress bar<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                block_list = list(red_ds.block_windows(1))<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                print(f\"Total blocks to process: {len(block_list)}\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                # Progress bar loop<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                for ji, window in tqdm(block_list, desc=\"Computing NDVI\", unit=\"block\"):<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    # Read both windows<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    red = red_ds.read(1, window=window).astype(\"float32\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    nir = nir_ds.read(1, window=window).astype(\"float32\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    # Nodata mask<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    mask = np.zeros(red.shape, dtype=bool)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    if red_nodata is not None:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                        mask |= (red == red_nodata)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    if nir_nodata is not None:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                        mask |= (nir == nir_nodata)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    denom = nir + red<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    safe = denom != 0<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    valid = safe &amp; (~mask)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    ndvi_block = np.full(red.shape, NDVI_NODATA, dtype=\"float32\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    ndvi_block[valid] = (nir[valid] - red[valid]) / denom[valid]<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    # Clamp<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    np.clip(ndvi_block, -1.0, 1.0, out=ndvi_block)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    out_ds.write(ndvi_block, 1, window=window)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>    total = (time.time() - t0) / 60<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>    print(f\"\\nFinished! Total time: {total:.2f} minutes\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>    print(f\"NDVI saved to: {OUTPUT_NDVI}\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>if __name__ == \"__main__\":<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>    compute_ndvi_streaming()<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>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 &lt;.1 were extracted and assigned the value 1. This is the trees layer.  The formula is listed below.<\/SPAN><\/P><P><SPAN /><\/P><P STYLE=\"text-indent:20;margin:1 1 1 20;\"><SPAN>SetNull(\"Mosaic_Height_Float2\", SetNull(\"NDVI_MERGED_Band_1\",1,\"Value &lt;.1\"),\"Value &lt; 3.048\")<\/SPAN><\/P><P><SPAN /><\/P><P STYLE=\"text-indent:20;\"><SPAN>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.<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>*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.<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN /><\/P><P STYLE=\"text-indent:20;\"><SPAN>See a few resources of research and municipalities using the 10ft / 3meter threshold:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><A href=\"https://urldefense.com/v3/__https://storymaps.arcgis.com/stories/51e0259c8daa42d6923e6a690754bfd6__;!!LVq-52eSbsHOEHc!8SwjW5BUipA0AJmisTVuNy78kc_B-Gm0v0xASgEWcktLvL96zPfNj3k9lP7I3UPJiJi2XV6HFol32-YUs_81eYvHWSuS1XcHeHxt$\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>Lake Oswego, OR<\/SPAN><\/A><\/P><P STYLE=\"text-indent:20;\"><A href=\"https://urldefense.com/v3/__https://www.cityservices.act.gov.au/trees-and-nature/trees/canopy-cover__;!!LVq-52eSbsHOEHc!8SwjW5BUipA0AJmisTVuNy78kc_B-Gm0v0xASgEWcktLvL96zPfNj3k9lP7I3UPJiJi2XV6HFol32-YUs_81eYvHWSuS1VVjRUeY$\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>Canberra, Australia<\/SPAN><\/A><\/P><P STYLE=\"text-indent:20;\"><A href=\"https://urldefense.com/v3/__https://vancouver.ca/files/cov/vancouver-tree-canopy-assessment-2022.pdf__;!!LVq-52eSbsHOEHc!8SwjW5BUipA0AJmisTVuNy78kc_B-Gm0v0xASgEWcktLvL96zPfNj3k9lP7I3UPJiJi2XV6HFol32-YUs_81eYvHWSuS1S4o_9JT$\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>Vancouver, Canada<\/SPAN><\/A><\/P><P STYLE=\"text-indent:20;\"><A href=\"https://urldefense.com/v3/__https://www.sciencedirect.com/science/article/pii/S2667393225000018__;!!LVq-52eSbsHOEHc!8SwjW5BUipA0AJmisTVuNy78kc_B-Gm0v0xASgEWcktLvL96zPfNj3k9lP7I3UPJiJi2XV6HFol32-YUs_81eYvHWSuS1Xtrdi-1$\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>https://www.sciencedirect.com/science/article/pii/S2667393225000018 <\/SPAN><\/A><SPAN>(section 2.3.1, includes a few references)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><A href=\"https://urldefense.com/v3/__https://www.mdpi.com/2072-4292/13/4/767__;!!LVq-52eSbsHOEHc!8SwjW5BUipA0AJmisTVuNy78kc_B-Gm0v0xASgEWcktLvL96zPfNj3k9lP7I3UPJiJi2XV6HFol32-YUs_81eYvHWSuS1QMlpC04$\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>https://www.mdpi.com/2072-4292/13/4/767 <\/SPAN><\/A><SPAN>(section 2.3.2)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>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.<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><A href=\"https://urldefense.com/v3/__https://treepeople.org/project/la-treecanopy-data/__;!!LVq-52eSbsHOEHc!8SwjW5BUipA0AJmisTVuNy78kc_B-Gm0v0xASgEWcktLvL96zPfNj3k9lP7I3UPJiJi2XV6HFol32-YUs_81eYvHWSuS1Y3qiHJP$\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>https://treepeople.org/project/la-treecanopy-data/<\/SPAN><\/A><\/P><P STYLE=\"text-indent:20;\"><A href=\"https://urldefense.com/v3/__https://storymaps.arcgis.com/stories/df083f2adb6a4650a738dbf2805674e2__;!!LVq-52eSbsHOEHc!8SwjW5BUipA0AJmisTVuNy78kc_B-Gm0v0xASgEWcktLvL96zPfNj3k9lP7I3UPJiJi2XV6HFol32-YUs_81eYvHWSuS1X2G4aK2$\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>https://storymaps.arcgis.com/stories/df083f2adb6a4650a738dbf2805674e2<\/SPAN><\/A><\/P><P STYLE=\"text-indent:20;\"><SPAN /><\/P><\/DIV><\/DIV><\/DIV>",
 "name": "LARIAC7/TREE_CANOPY_2024",
 "description": "<DIV STYLE=\"text-align:Left;\"><DIV><DIV><P><SPAN>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. <\/SPAN><\/P><P><SPAN /><\/P><P STYLE=\"text-indent:20;\"><SPAN>1. <\/SPAN><A href=\"https://lariac-lacounty.hub.arcgis.com/\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>Los Angeles Region Imagery Acquisition Consortium (LARIAC-7)  <\/SPAN><\/A><SPAN>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. <\/SPAN><A href=\"https://gisgeography.com/ndvi-normalized-difference-vegetation-index/\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>Click here for more information.<\/SPAN><\/A><\/P><P><SPAN /><\/P><P STYLE=\"text-indent:20;\"><SPAN>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.<\/SPAN><\/P><P><SPAN /><\/P><P><SPAN STYLE=\"font-weight:bold;\">Processing Steps:<\/SPAN><\/P><P><SPAN>1. To create NDVI, the following code was used outside of ArcPro to efficiently combine the Red and Infrared bands of the CIR imagery.<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>import os<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>import time<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>import numpy as np<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>import rasterio<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>from tqdm import tqdm   # progress bar<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN># =========================<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN># CONFIG SECTION (EDIT ME)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN># =========================<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>RED_RASTER = r\"Red.tif\"<\/SPAN><\/P><P \n\t\t\t\tSTYLE=\"text-indent:20;\"><SPAN>NIR_RASTER = r\"Infrared.tif\"<\/SPAN><\/P><P \n\t\t\t\tSTYLE=\"text-indent:20;\"><SPAN>OUTPUT_NDVI = r\"NDVI_2024_V2.tif\"<\/SPAN><\/P><P \n\t\t\t\tSTYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>NDVI_NODATA = -9999.0<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>OUTPUT_COMPRESS = \"NONE\"        # FASTEST. Switch to DEFLATE later if needed.<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>BLOCK_SIZE = 512                # Good block alignment for most rasters.<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>os.environ.setdefault(\"GDAL_NUM_THREADS\", \"ALL_CPUS\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>def compute_ndvi_streaming():<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>    t0 = time.time()<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>    print(f\"\\nOpening RED raster: {RED_RASTER}\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>    with rasterio.open(RED_RASTER) as red_ds:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>        w = red_ds.width<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>        h = red_ds.height<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>        red_nodata = red_ds.nodata<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>        print(f\"  Size: {w} × {h}\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>        <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>        with rasterio.open(NIR_RASTER) as nir_ds:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            if (nir_ds.width, nir_ds.height) != (w, h):<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                raise ValueError(\"Input rasters differ in size.\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            if nir_ds.transform != red_ds.transform:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                raise ValueError(\"Input rasters have different transforms.\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            if nir_ds.crs != red_ds.crs:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                raise ValueError(\"Input rasters differ in CRS.\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            nir_nodata = nir_ds.nodata<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            profile = red_ds.profile.copy()<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            profile.update(<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                driver=\"GTiff\",<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                dtype=\"float32\",<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                count=1,<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                nodata=NDVI_NODATA,<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                tiled=True,<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                blockxsize=BLOCK_SIZE,<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                blockysize=BLOCK_SIZE,<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                compress=OUTPUT_COMPRESS,<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                BIGTIFF=\"IF_SAFER\",<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            )<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            print(\"\\nCreating output NDVI raster...\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>            with rasterio.open(OUTPUT_NDVI, \"w\", **profile) as out_ds:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                # Count blocks for progress bar<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                block_list = list(red_ds.block_windows(1))<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                print(f\"Total blocks to process: {len(block_list)}\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                # Progress bar loop<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                for ji, window in tqdm(block_list, desc=\"Computing NDVI\", unit=\"block\"):<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    # Read both windows<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    red = red_ds.read(1, window=window).astype(\"float32\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    nir = nir_ds.read(1, window=window).astype(\"float32\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    # Nodata mask<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    mask = np.zeros(red.shape, dtype=bool)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    if red_nodata is not None:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                        mask |= (red == red_nodata)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    if nir_nodata is not None:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                        mask |= (nir == nir_nodata)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    denom = nir + red<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    safe = denom != 0<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    valid = safe &amp; (~mask)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    ndvi_block = np.full(red.shape, NDVI_NODATA, dtype=\"float32\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    ndvi_block[valid] = (nir[valid] - red[valid]) / denom[valid]<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    # Clamp<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    np.clip(ndvi_block, -1.0, 1.0, out=ndvi_block)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>                    out_ds.write(ndvi_block, 1, window=window)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>    total = (time.time() - t0) / 60<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>    print(f\"\\nFinished! Total time: {total:.2f} minutes\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>    print(f\"NDVI saved to: {OUTPUT_NDVI}\")<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN> <\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>if __name__ == \"__main__\":<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>    compute_ndvi_streaming()<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>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 &lt;.1 were extracted and assigned the value 1. This is the trees layer.  The formula is listed below.<\/SPAN><\/P><P><SPAN /><\/P><P STYLE=\"text-indent:20;margin:1 1 1 20;\"><SPAN>SetNull(\"Mosaic_Height_Float2\", SetNull(\"NDVI_MERGED_Band_1\",1,\"Value &lt;.1\"),\"Value &lt; 3.048\")<\/SPAN><\/P><P><SPAN /><\/P><P STYLE=\"text-indent:20;\"><SPAN>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.<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>*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.<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN /><\/P><P STYLE=\"text-indent:20;\"><SPAN>See a few resources of research and municipalities using the 10ft / 3meter threshold:<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><A href=\"https://urldefense.com/v3/__https://storymaps.arcgis.com/stories/51e0259c8daa42d6923e6a690754bfd6__;!!LVq-52eSbsHOEHc!8SwjW5BUipA0AJmisTVuNy78kc_B-Gm0v0xASgEWcktLvL96zPfNj3k9lP7I3UPJiJi2XV6HFol32-YUs_81eYvHWSuS1XcHeHxt$\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>Lake Oswego, OR<\/SPAN><\/A><\/P><P STYLE=\"text-indent:20;\"><A href=\"https://urldefense.com/v3/__https://www.cityservices.act.gov.au/trees-and-nature/trees/canopy-cover__;!!LVq-52eSbsHOEHc!8SwjW5BUipA0AJmisTVuNy78kc_B-Gm0v0xASgEWcktLvL96zPfNj3k9lP7I3UPJiJi2XV6HFol32-YUs_81eYvHWSuS1VVjRUeY$\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>Canberra, Australia<\/SPAN><\/A><\/P><P STYLE=\"text-indent:20;\"><A href=\"https://urldefense.com/v3/__https://vancouver.ca/files/cov/vancouver-tree-canopy-assessment-2022.pdf__;!!LVq-52eSbsHOEHc!8SwjW5BUipA0AJmisTVuNy78kc_B-Gm0v0xASgEWcktLvL96zPfNj3k9lP7I3UPJiJi2XV6HFol32-YUs_81eYvHWSuS1S4o_9JT$\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>Vancouver, Canada<\/SPAN><\/A><\/P><P STYLE=\"text-indent:20;\"><A href=\"https://urldefense.com/v3/__https://www.sciencedirect.com/science/article/pii/S2667393225000018__;!!LVq-52eSbsHOEHc!8SwjW5BUipA0AJmisTVuNy78kc_B-Gm0v0xASgEWcktLvL96zPfNj3k9lP7I3UPJiJi2XV6HFol32-YUs_81eYvHWSuS1Xtrdi-1$\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>https://www.sciencedirect.com/science/article/pii/S2667393225000018 <\/SPAN><\/A><SPAN>(section 2.3.1, includes a few references)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><A href=\"https://urldefense.com/v3/__https://www.mdpi.com/2072-4292/13/4/767__;!!LVq-52eSbsHOEHc!8SwjW5BUipA0AJmisTVuNy78kc_B-Gm0v0xASgEWcktLvL96zPfNj3k9lP7I3UPJiJi2XV6HFol32-YUs_81eYvHWSuS1QMlpC04$\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>https://www.mdpi.com/2072-4292/13/4/767 <\/SPAN><\/A><SPAN>(section 2.3.2)<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><SPAN>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.<\/SPAN><\/P><P STYLE=\"text-indent:20;\"><A href=\"https://urldefense.com/v3/__https://treepeople.org/project/la-treecanopy-data/__;!!LVq-52eSbsHOEHc!8SwjW5BUipA0AJmisTVuNy78kc_B-Gm0v0xASgEWcktLvL96zPfNj3k9lP7I3UPJiJi2XV6HFol32-YUs_81eYvHWSuS1Y3qiHJP$\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>https://treepeople.org/project/la-treecanopy-data/<\/SPAN><\/A><\/P><P STYLE=\"text-indent:20;\"><A href=\"https://urldefense.com/v3/__https://storymaps.arcgis.com/stories/df083f2adb6a4650a738dbf2805674e2__;!!LVq-52eSbsHOEHc!8SwjW5BUipA0AJmisTVuNy78kc_B-Gm0v0xASgEWcktLvL96zPfNj3k9lP7I3UPJiJi2XV6HFol32-YUs_81eYvHWSuS1X2G4aK2$\" target=\"_blank\" STYLE=\"text-decoration:underline;\"><SPAN>https://storymaps.arcgis.com/stories/df083f2adb6a4650a738dbf2805674e2<\/SPAN><\/A><\/P><P STYLE=\"text-indent:20;\"><SPAN /><\/P><\/DIV><\/DIV><\/DIV>",
 "extent": {
  "xmin": 6273759.205987,
  "ymin": 1564795.154747,
  "xmax": 6672403.381318,
  "ymax": 2100719.2252970003,
  "spatialReference": {
   "wkid": 103007,
   "latestWkid": 6424
  }
 },
 "initialExtent": {
  "xmin": 6273759.205987,
  "ymin": 1564795.154747,
  "xmax": 6672403.381318,
  "ymax": 2100719.2252970003,
  "spatialReference": {
   "wkid": 103007,
   "latestWkid": 6424
  }
 },
 "fullExtent": {
  "xmin": 6273759.205987,
  "ymin": 1564795.154747,
  "xmax": 6672403.381318,
  "ymax": 2100719.2252970003,
  "spatialReference": {
   "wkid": 103007,
   "latestWkid": 6424
  }
 },
 "hasMultidimensions": false,
 "pixelSizeX": 3.280833000000002,
 "pixelSizeY": 3.2808330000000017,
 "datasetFormat": "FGDBR",
 "uncompressedSize": 19848168450,
 "blockWidth": 128,
 "blockHeight": 128,
 "compressionType": "LZ77",
 "bandNames": [
  "Band_1"
 ],
 "allowCopy": true,
 "allowAnalysis": true,
 "bandCount": 1,
 "pixelType": "U8",
 "minPixelSize": 0,
 "maxPixelSize": 0,
 "copyrightText": "Public Works SMPM, Nathan Shapiro",
 "serviceDataType": "esriImageServiceDataTypeGeneric",
 "serviceSourceType": "esriImageServiceSourceTypeDataset",
 "minValues": [
  1
 ],
 "maxValues": [
  1
 ],
 "meanValues": [
  1
 ],
 "stdvValues": [
  0
 ],
 "objectIdField": "",
 "fields": [],
 "capabilities": "Image,Metadata,Mensuration",
 "defaultMosaicMethod": "Center",
 "allowedMosaicMethods": "",
 "sortField": "",
 "sortValue": null,
 "sortAscending": true,
 "mosaicOperator": "First",
 "defaultCompressionQuality": 75,
 "defaultResamplingMethod": "Bilinear",
 "maxImageHeight": 4100,
 "maxImageWidth": 15000,
 "allowRasterFunction": true,
 "rasterFunctionInfos": [
  {
   "name": "None",
   "description": "",
   "help": ""
  }
 ],
 "rasterTypeInfos": [
  {
   "name": "Raster Dataset",
   "description": "Supports all ArcGIS Raster Datasets",
   "help": ""
  }
 ],
 "mensurationCapabilities": "Basic",
 "hasHistograms": true,
 "hasColormap": false,
 "hasRasterAttributeTable": true,
 "minScale": 0,
 "maxScale": 0,
 "exportTilesAllowed": false,
 "supportsStatistics": false,
 "supportsAdvancedQueries": false,
 "editFieldsInfo": null,
 "ownershipBasedAccessControlForRasters": null,
 "allowComputeTiePoints": false,
 "useStandardizedQueries": true,
 "advancedQueryCapabilities": {
  "useStandardizedQueries": true,
  "supportsStatistics": false,
  "supportsOrderBy": false,
  "supportsDistinct": false,
  "supportsPagination": false
 },
 "spatialReference": {
  "wkid": 103007,
  "latestWkid": 6424
 }
}