R

vegsub <- veg2017[which(veg2017$adjEasting >= ext(chm)[1] &
                        veg2017$adjEasting <= ext(chm)[2] &
                        veg2017$adjNorthing >= ext(chm)[3] & 
                        veg2017$adjNorthing <= ext(chm)[4]),]

Python

vegsub = veg2017.loc[(veg2017["adjEasting"] >= chm.bounds[0]) &
                 (veg2017["adjEasting"] <= chm.bounds[1]) &
                 (veg2017["adjNorthing"] >= chm.bounds[2]) & 
                 (veg2017["adjNorthing"] <= chm.bounds[3])]
vegsub = vegsub.reset_index(drop=True)

R

The extract() function from the terra package gets the values from the tile at the given coordinates.

valCHM <- extract(chm, 
                  cbind(vegsub$adjEasting,
                  vegsub$adjNorthing))
plot(valCHM$NEON_D16_WREF_DP3_580000_5075000_CHM~
vegsub$height, pch=20, xlab="Height",
ylab="Canopy height model")
lines(c(0,50), c(0,50), col="grey")

How strong is the correlation between the ground and lidar measurements?

cor(valCHM$NEON_D16_WREF_DP3_580000_5075000_CHM, 
    vegsub$height, use="complete")

## [1] 0.4036681

Python

The sample_gen() function from the rasterio.sample module gets the values from the tile at the given coordinates.

valCHM = list(sample.sample_gen(chm, 
                                tuple(zip(vegsub["adjEasting"], 
                                          vegsub["adjNorthing"])),
                                masked=True))
fig, ax = plt.subplots()
ax.plot((0,50), (0,50), linewidth=1, color="black")
ax.scatter(vegsub.height, valCHM, s=1)
ax.set_xlabel("Height")
ax.set_ylabel("Canopy height model")
plt.show()

How strong is the correlation between the ground and lidar measurements?

CHMlist = np.array([c.tolist()[0] for c in valCHM])
idx = np.intersect1d(np.where(np.isfinite(vegsub.height)), 
                     np.where(CHMlist != None))
np.corrcoef(vegsub.height[idx], list(CHMlist[idx]))[0,1]

## 0.40366813643317934

R

valCHMbuff <- extract(chm, 
                  buffer(vect(cbind(vegsub$adjEasting,
                  vegsub$adjNorthing)),
                  width=vegsub$adjCoordinateUncertainty),
                  fun=max)
plot(valCHMbuff$NEON_D16_WREF_DP3_580000_5075000_CHM~
vegsub$height, pch=20, xlab="Height",
ylab="Canopy height model")
lines(c(0,50), c(0,50), col="grey")

cor(valCHMbuff$NEON_D16_WREF_DP3_580000_5075000_CHM, 
    vegsub$height, use="complete")

## [1] 0.3903361

Python

To extract values using a buffer in Python, we need to create a shapefile of the buffered locations, and then extract the maximum value for each area in the shapefile.

vegloc = vegsub[["individualID","adjEasting","adjNorthing","adjCoordinateUncertainty"]]
v = vegloc.rename(columns={"individualID": "indID", "adjEasting": "easting",
               "adjNorthing": "northing", "adjCoordinateUncertainty": "coordUnc"},
               inplace=False)
gdf = gpd.GeoDataFrame(
v, geometry=gpd.points_from_xy(v.easting, v.northing))
gdf["geometry"] = gdf["geometry"].buffer(distance=gdf["coordUnc"])
gdf.to_file(os.getcwd() + "/trees_with_buffer.shp")
chm_height = rs.zonal_stats(os.getcwd() + "/trees_with_buffer.shp", chmx.values,
affine=chmx.rio.transform(),
nodata=-9999, stats="max")
valCHMbuff = [h["max"] for h in chm_height]

And plot the results:

fig, ax = plt.subplots()
ax.plot((0,50), (0,50), linewidth=1, color="black")
ax.scatter(vegsub.height, valCHMbuff, s=1)
ax.set_xlabel("Height")
ax.set_ylabel("Canopy height model")
plt.show()

CHMbufflist = np.array(valCHMbuff)
idx = np.intersect1d(np.where(np.isfinite(vegsub.height)), 
                     np.where(CHMbufflist != None))
np.corrcoef(vegsub.height[idx], list(CHMbufflist[idx]))[0,1]

## 0.39033606753955524

R

easting10 <- 10*floor(vegsub$adjEasting/10)
northing10 <- 10*floor(vegsub$adjNorthing/10)
vegsub <- cbind(vegsub, easting10, northing10)

Use the aggregate() function to get the tallest tree in each 10m bin.

vegbin <- stats::aggregate(vegsub, 
                           by=list(vegsub$easting10, 
                                   vegsub$northing10), 
                           FUN=max)

To get the CHM values for the 10m bins, use the terra package version of the aggregate() function. Let’s take a look at the lower-resolution image we get as a result.

CHM10 <- terra::aggregate(chm, fact=10, fun=max)
plot(CHM10, col=topo.colors(5))

Use the extract() function again to get the values from each pixel. Our grids are numbered by the corners, so add 5 to each tree coordinate to make sure it’s in the correct pixel.

vegbin$easting10 <- vegbin$easting10 + 5
vegbin$northing10 <- vegbin$northing10 + 5
binCHM <- extract(CHM10, cbind(vegbin$easting10, 
                               vegbin$northing10))
plot(binCHM$NEON_D16_WREF_DP3_580000_5075000_CHM~
       vegbin$height, pch=20, 
     xlab="Height", ylab="Canopy height model")
lines(c(0,50), c(0,50), col="grey")

cor(binCHM$NEON_D16_WREF_DP3_580000_5075000_CHM, 
    vegbin$height, use="complete")

## [1] 0.3472434

Python

vegsub["easting10"] = 10*np.floor(vegsub.adjEasting/10)
vegsub["northing10"] = 10*np.floor(vegsub.adjNorthing/10)
vegsubloc = vegsub[["height","easting10","northing10"]]

Use the groupby() function to get the tallest tree in each 10m bin.

vegbin = vegsubloc.groupby(["easting10", "northing10"]).max().add_suffix('_max').reset_index()

To get the CHM values for the 10m bins, use the warp.reproject() method in the rasterio package.

target_res = (10, 10)
with rasterio.open(os.getcwd() + "/DP3.30015.001/neon-aop-products/2017/FullSite/D16/2017_WREF_1/L3/DiscreteLidar/CanopyHeightModelGtif/NEON_D16_WREF_DP3_580000_5075000_CHM.tif") as src:
data, transform = rasterio.warp.reproject(source=src.read(),
src_transform=src.transform,
src_crs=src.crs,
dst_crs=src.crs,
dst_nodata=src.nodata,
dst_resolution=target_res,
resampling=Resampling.max)
profile = src.profile
profile.update(transform=transform, driver='GTiff',
height=data.shape[1], width=data.shape[2])
 with rasterio.open(os.getcwd() + '/CHM_10m.tif', 'w', **profile) as dst:
                dst.write(data)

chm10 = rasterio.open(os.getcwd() + '/CHM_10m.tif')

Let’s take a look at the lower-resolution image we get as a result.

plt.imshow(chm10.read(1))
plt.show()

Use the sample() function again to get the values from the pixel corresponding to each maximum tree height estimate. Our grids are numbered by the corners, so add 5 to each tree coordinate to make sure it’s in the correct pixel.

valCHM10 = list(sample.sample_gen(chm10, tuple(zip(vegbin["easting10"]+5,
                                                   vegbin["northing10"]+5)),
                                                   masked=True))
fig, ax = plt.subplots()
ax.plot((0,50), (0,50), linewidth=1, color="black")
ax.scatter(vegbin.height_max, valCHM10, s=1)
ax.set_xlabel("Height")
ax.set_ylabel("Canopy height model")
plt.show()

CHM10list = np.array([c.tolist()[0] for c in valCHM10])
idx = np.intersect1d(np.where(np.isfinite(vegbin.height_max)), 
                     np.where(CHM10list != None))
np.corrcoef(vegbin.height_max[idx], list(CHM10list[idx]))[0,1]

## 0.346718246893291

R

vegsub <- vegsub[order(vegsub$height, 
                       decreasing=T),]

Now, for each tree, let’s estimate which nearby trees might be beneath its canopy, and discard those points. To do this:

1. Calculate the distance of each tree from the target tree.
2. Pick a reasonable estimate for canopy size, and discard shorter trees within that radius. The radius I used is 0.3 times the height, based on some rudimentary googling about Douglas fir allometry. It could definitely be improved on!
3. Iterate over all trees.

We’ll use a simple for loop to do this:

vegfil <- vegsub
for(i in 1:nrow(vegsub)) {
    if(is.na(vegfil$height[i]))
        next
    dist <- sqrt((vegsub$adjEasting[i]-vegsub$adjEasting)^2 + 
                (vegsub$adjNorthing[i]-vegsub$adjNorthing)^2)
    vegfil$height[which(dist<0.3*vegsub$height[i] & 
                        vegsub$height<vegsub$height[i])] <- NA
}
vegfil <- vegfil[which(!is.na(vegfil$height)),]

Now extract the raster values, as above.

filterCHM <- extract(chm, 
                     cbind(vegfil$adjEasting, 
                           vegfil$adjNorthing))
plot(filterCHM$NEON_D16_WREF_DP3_580000_5075000_CHM~
       vegfil$height, pch=20, 
     xlab="Height", ylab="Canopy height model")
lines(c(0,50), c(0,50), col="grey")

cor(filterCHM$NEON_D16_WREF_DP3_580000_5075000_CHM,
    vegfil$height)

## [1] 0.8377803

Python

vegfil = vegsub.sort_values(by="height", ascending=False, ignore_index=True)

Now, for each tree, let’s estimate which nearby trees might be beneath its canopy, and discard those points. To do this:

1. Calculate the distance of each tree from the target tree.
2. Pick a reasonable estimate for canopy size, and discard shorter trees within that radius. The radius I used is 0.3 times the height, based on some rudimentary googling about Douglas fir allometry. It could definitely be improved on!
3. Iterate over all trees.

We’ll use a simple for loop to do this:

height = vegfil.height.reset_index()
for i in vegfil.index:
    if height.height[i] is None:
        pass
    else:
        dist = np.sqrt(np.square(vegfil.adjEasting[i]-vegfil.adjEasting) + 
                       np.square(vegfil.adjNorthing[i]-vegfil.adjNorthing))
        idx = vegfil.index[(vegfil.height<height.height[i]) & (dist<0.3*height.height[i])]
        height.loc[idx, "height"] = None

Now extract the raster values, as above.

filterCHM = list(sample.sample_gen(chm, tuple(zip(vegfil["adjEasting"],
                                                  vegfil["adjNorthing"])),
                                                  masked=True))
fig, ax = plt.subplots()
ax.plot((0,50), (0,50), linewidth=1, color="black")
ax.scatter(height.height, filterCHM, s=1)
ax.set_xlabel("Height")
ax.set_ylabel("Canopy height model")
plt.show()

filCHMlist = np.array([c.tolist()[0] for c in filterCHM])
idx = np.intersect1d(np.where(np.isfinite(height.height)), 
                     np.where(filCHMlist != None))
np.corrcoef(height.height[idx], list(filCHMlist[idx]))[0,1]

## 0.8380927468867392

R

vegfil <- vegfil[which(vegfil$plantStatus=="Live"),]
filterCHM <- extract(chm,
cbind(vegfil$adjEasting,
vegfil$adjNorthing))
plot(filterCHM$NEON_D16_WREF_DP3_580000_5075000_CHM~
vegfil$height, pch=20,
xlab="Height", ylab="Canopy height model")
lines(c(0,50), c(0,50), col="grey")

cor(filterCHM$NEON_D16_WREF_DP3_580000_5075000_CHM,
    vegfil$height)

## [1] 0.9145612

Python

idx = vegfil.index[vegfil.plantStatus!="Live"]
height.loc[idx, "height"] = None
fig, ax = plt.subplots()
ax.plot((0,50), (0,50), linewidth=1, color="black")
ax.scatter(height.height, filterCHM, s=1)
ax.set_xlabel("Height")
ax.set_ylabel("Canopy height model")
plt.show()

idx = np.intersect1d(np.where(np.isfinite(height.height)), 
                     np.where(filCHMlist != None))
np.corrcoef(height.height[idx], list(filCHMlist[idx]))[0,1]

## 0.9146181359857672