PatchMorph Issue—Gap and Spur Thresholds, 0s and 1s Mis-Assigned?

[paigekouba]

Hi there! I am hoping to use the patchMorph function to delineate openings in a set of forest plots. I am having two problems:

1) the "suitable" (1) and "unsuitable" (0) values seem to be getting swapped somewhere along the way; 2) the gap and spur threshold values I've set don't lead to reasonable results.

I start with a spatRaster of my 1-ha plot, with all crown areas assigned "0" (unsuitable) and all other areas assigned "1" (suitable), and crs set to "local" (ie a Cartesian plane with units in m—our data aren't mapped with GIS, just X and Y coords). I tried lots of combinations of gap and spur; one that produced reasonable-looking opening polygons was 8 and 4. However, I understand that the gap threshold should be about the size of the smallest crowns in my dataset, and the spur threshold should be about the size of the largest crowns in my dataset. If I used 2 and 10, which is more consistent with my dataset, the ecological meanings of the gap and spur thresholds, and other published studies, I get a plot with almost no crown area.

Here is the most relevant excerpt of my code: plot(sr_test, col=c("white", "yellow"), main="Starting SpatRaster") # 0s in crowns (white), 1s everywhere else (yellow) pm.rast <- patchMorph(sr_test, buffer = 5, suitThresh = 1, gapThresh = 8, spurThresh = 4, verbose = TRUE) plot(pm.rast, main="PatchMorph Results (Gap-10 & Spur-6)", col=c("yellow", "forestgreen")) # now the areas near crowns are 1 and openings are 0! ? plot(crowns, col="black", add=TRUE)

Do you know what I might be doing wrong? I have considered that there might be an issue with converting the pixels (in the circular kernel) to meters (to match my data); alternatively I might be getting my 0s and 1s confused; or I may just be misunderstanding the parameters themselves. So far, the threshold values seem just plain funky. I would greatly appreciate any advice on how to get them to behave! Reproducible code attached. PatchMorph_Issue_pvk.txt

Thank you, Paige

P.S. Interestingly, for just one of my 12 plots, patchMorph doesn't behave as described above. In that case, the shapes drawn by the algorithm seem to be focused on the crown areas, not the gaps; however if you ignore that rather significant detail, the 2/10 gap/spur threshold settings produce pretty decent shapes in this one case! I'm totally flummoxed.

[bi0m3trics]

Hi @paigekouba, I appreciate the issue, and I'll get back on this as soon as I can. I have a ton of traveling this month, so I'll get to this/back on this ASAP. Some of this seems like I need to document the package and functions better (no one's used it!). Couple questions:

1. I can confirm the switching is happening, and I'm not 100% sure why yet. The threshold argument is so habitat suitability can be specified for a range of raster values (e.g., 0-3 is unsuitable, and 3-27 is suitable) and I'm not really sure of how to make it work for a "standardized" case. For example, if threshold is 1, then I can assume above 1 is suitable and below 1 is unsuitable, however, when someone reverses the classification and the threshold is still 1, I'm trying to figure out how to handle that case. Maybe I just need to think about it more....
2. When you say, "threshold values I've set don't lead to reasonable results," is this based on gut, hand calculations, or comparisons to some other package or code? Again, I need to document things better, but the gap and spurs values are in resolution of the raster, i.e., your raster is 0.1 m (based on code provided) so a spur of 10 would be 1m in your raster's case.

I may have some time while traveling later this week, so I'll provide an update if I find a solution or confirm the code is running as intended. Andrew

[paigekouba]

Thanks very much for your help with this, especially while you're on the road!

1. I see how the binary openings-vs-crown case doesn't translate well to the suitability threshold. Maybe I can try setting an arbitrary higher value for the "suitable" areas, and see what happens?
2. Based on other papers I've read using the ArcGIS version of PatchMorph (J. Lydersen 2013, D. Fry 2014, J. Ng 2020), the gap and spur values should be set relative to the min and max crown diameters in the dataset, respectively, so I had been expecting I would use values around 2m for gap and 10m for spur. Even if I correct for the resolution, and use e.g. 20 and 100 for those parameters, it takes a very long time to run and the polygons it draws don't fit the scale of the crowns I'm trying to map openings for. Attaching a photo of this to demonstrate

[bi0m3trics]

@paigekouba Ok, so I had some time to get back into this (I was traveling quite a bit), and here's my current takeaways:

1. The first parameter as defined in the PatchMorth paper, is a density threshold and thus defines the minimum density of suitable habitat required for an area to be considered as part of a habitat patch. It is used in the initial step of the PatchMorph algorithm, where a suitable land cover neighborhood density filter excludes areas of low density suitable habitat from being part of a patch. This threshold is important because it helps ensure that the delineated patches represent areas with a significant presence of suitable habitat, reflecting the habitat requirements of the focus organism. Additionally, it makes sense that the provided raster is either binary (0,1) or increasing numerical, (0,1,2,3...). The first image below is an example using gap and spur thresholds set to 6 overlain with your tree crowns. This look right to me...
2. I don't have access to what others have used to calculate PathcMorph rasters (Girvetz and Greco's 2007 ArcGIS extension), so I can't comment to their accuracy. That's why I created this package, and it's the strictest interpretation of their algorithm based on how I read the paper. If someone has access to the original code, I'd be happy to try and convert and compare. I've seen script others have used, but I wasn't convinced the concessions they made were correct (sums instead of maximums in the focal functions).
3. It looks like terra (which I didn't know) does everything using the units for the crs, and makes any conversions necessary behind the scenes. Thus, if you have a crs or you specify one, the spur and gap parameters will be in those units (which have to be the same). I never noticed this when I converted form raster to terra
4. I added/revised the help and documentation, so that it may be clearer to the user what the parameters should be...let me know if things are still not clear and I'll revise further.
5. Makes sure to be careful when plotting. I think one of the earlier images you showed me, was the result one class (0 or 1) not being present at huge gaps and spurs, and thus the coloring of the palette shifted down.
6. Everything seems to be working as intended.

PatchMorph Results (Gap-6m & Spur-6m)

PatchMorph Results for Gap& Spur seq(1, 12, by=2) or c(1 ,3 ,5 ,7 ,9,11)

-------------------------------------------

library(sf)
library(terra)

# X values
Xs <- c(-54.3, -53.6, -52.4, -48.5, -47.0, -45.8, -45.6, -40.6, -39.5, -39.2,
        -36.1, -35.2, -34.6, -34.1, -30.3, -28.2, -27.7, -26.0, -25.9, -25.9,
        -23.1, -21.6, -21.2, -20.5, -17.9, -17.3, -15.9, -15.4, -14.6, -10.0,
        -5.5, -4.9, -2.3, -1.9, 0.3, 2.2, 4.3, 4.7, 6.0, 6.9, 8.6, 9.1, 10.4,
        13.5, 13.9, 16.8, 22.3, 28.1, 31.4, 36.0, 37.7, 50.0, 50.7, 5.2, -3.3,
        27.4, 26.6, -2.2, -5.4, 0.6, -29.9)
# Y values
Ys <- c(3.2, -5.1, 7.2, -20.5, -25.3, -26.0, -21.3, -0.2, 12.7, -34.4, -35.6,
        36.5, 23.5, -9.3, -1.9, 35.4, 18.0, -37.6, -26.9, -11.9, -41.0, 24.3,
        12.2, 30.0, 5.4, 32.9, -4.3, -28.8, 24.1, -23.0, -55.7, 14.4, 12.2,
        -53.1, -7.1, -13.8, -13.2, 48.1, 6.5, 35.5, 10.4, 41.7, 16.7, 53.4,
        52.0, -12.8, -30.8, -25.0, -40.9, -22.7, 41.5, -21.3, 12.3, 20.3,
        23.6, 40.1, -10.8, -32.3, -34.9, -36.0, 24.0)
# Crown values
Cs <- c(1.4116, 1.1624, 1.2336, 3.5476, 2.4084, 1.7676, 2.7644, 1.1268, 3.7612,
        1.5540, 1.2336, 3.3340, 1.1980, 1.1624, 2.2304, 2.6220, 3.6188, 1.5540,
        3.4764, 3.8324, 1.5896, 2.4440, 1.1268, 3.2272, 3.7968, 3.4764, 3.5476,
        3.7256, 4.1528, 2.7644, 1.5896, 2.8000, 4.2240, 1.9812, 4.4732, 3.4408,
        2.8356, 1.2336, 3.6188, 1.4828, 3.9392, 1.8744, 5.1140, 3.2984, 3.5120,
        3.9748, 3.4408, 3.0136, 4.1528, 5.2208, 4.0460, 2.6932, 4.3308, 2.9424,
        1.5896, 1.5540, 4.0104, 2.1236, 2.2304, 2.6576, 1.5540)
# Create data frame
pm_df <- data.frame(X = Xs, Y = Ys, crown = Cs)

ctr = data.frame(X = 0, Y = 0) # plot center to draw boundary of 1ha circle
bound = st_as_sf(ctr, coords = c("X", "Y")) |> st_buffer(sqrt(10000/pi)) # boundary of 1ha circle
stems <- st_as_sf(pm_df, coords = c("X", "Y")) # points for each tree w/dbh attribute
crowns = st_buffer(stems, dist = stems$crown) # each crown defined as a circle with r=crown radius
notcrowns <- st_difference(bound, st_union(crowns)) # all area in non-crown space
yescrowns <- st_difference(bound, st_union(notcrowns)) # all area in crown space

r_notcrowns <- rasterize(notcrowns, rast(ext(bound), res = 0.1)) # turns sf into a rasterlayer -- this one is for all non-crown areas
# want now to make a raster with 1s for not-crown and 0s for crown
r_yescrowns <- rasterize(yescrowns, rast(ext(bound), res = 0.1)) # sf --> rasterlayer, this time for all crown areas

values(r_yescrowns)[values(r_yescrowns) == 1] <- 0 # reassign all values in crown areas to 0 ("unsuitable")

sr_test <- merge(r_notcrowns, r_yescrowns) # creates one rasterlayer with 0s in crown areas, 1s everywhere else
crs(sr_test) <- "local" # set crs to Cartesian plane in meters

# flip it so that 1s are in crown areas and 0s are everywhere else
sr_test <- 1-sr_test

library(patchwoRk)

# single run
par(mfrow=c(1,2))
plot(sr_test, main="Starting SpatRaster")
plot(crowns, col=adjustcolor("black", alpha=0.3), add=TRUE)
pm.rast <- patchMorph(sr_test, buffer = 6, suitThresh = 1, gapThresh = 6, spurThresh = 6, verbose = TRUE)
plot(pm.rast,  main="PatchMorph Results (Gap-6m & Spur-6m)")
plot(crowns, col=adjustcolor("black", alpha=0.3), add=TRUE)

# multi-run
par(mfrow=c(6,6))
for (i in seq(1, 12, by=2)) {
 for (j in seq(1, 12, by=2)) {
   pm.rast <- patchMorph(sr_test, buffer = max(i,j)*2, suitThresh = 1, gapThresh = i, spurThresh = j, verbose = TRUE)
   plot(pm.rast,  main=paste("Gap-", i, "m & Spur-", j, "m"))
   plot(crowns, col=adjustcolor("black", alpha=0.3), add=TRUE)
 }
}

Give this a run and let me know if it looks like it's working. There's still some funky edge issues I'd like to resolve. but that'll have to wait a bit.

[paigekouba]

Hi Andrew, thanks again, your comments are helpful! I'll respond to your points in order:

1. Ok that makes sense about the suitability threshold, thank you for explaining
2. I actually know someone who just finished his PhD with Steve Greco, would you like me to connect you, or see if I can get a copy of the PatchMorph code from their 2007 paper?
3. Ok good on the units issue! My lab mate told me about terra, I understand it's fairly new but has a lot of good features. Having never used raster I can't compare, but it seems pretty useful to me
4. Documentation looks good, thanks!
5. Duly noted!
6. Can you comment on how you would choose gap and spur thresholds, a priori? My plots match yours now when I use G-6 S-6 (or any of the combinations in your for loop), so that's good. However, the three papers I know using PatchMorph (Jamie Lydersen 2013, Danny Fry 2014, and Jan Ng 2020) all used gap thresholds from 1.5-3m, based on the minimum tree crown diameter in their study; they all used spur thresholds from 10-12m, based on the maximum tree crown diameter. If I use the same logic, I would have a gapThresh of 2m and a spurThresh of 10m, but as shown in your plots those parameters yield almost 100% gap area. Right now I think the shapes drawn by G-9 S-7 are the most informative, and look the most like the shapes in those other papers, but I don't have a good (eco)logical explanation for choosing them.

If it would be easier to talk about this synchronously, perhaps we could discuss via Zoom? Let me know if you might have time for that, and I can email you! Thanks, Paige

[bi0m3trics]

@paigekouba I'd start with min and max crown radii... or the radius of expected individual tree cover/openings (assuming circular areas). If Danny and Jamie used 3,12; then try starting with 1.5, 6. -a

[paigekouba]

Hi again Andrew, I just noticed something that might interest you! For my dataset, the minimum crown radius is 0.74m (10th percentile is 1.25). Max radius is 6.0m (90th percentile is 5.0). Based on those, I would expect to have a gapThresh of ~2m and a spurThresh of ~10m. But as you can see from the 6x6 grid of maps you provided with the different parameter combos, 2 and 10 (or 1 and 9, which round up to the same) gives almost 100% openings: in fact, that map has only got a few green spots inside of areas that should be occupied by tree crowns. Meanwhile, the map for gap-11 and spur-9 draws openings (gray blobs) that look reasonable for the tree crowns in the sample plot. Certain other combos (9/7, 7/5, 5/3) also look pretty good, just with a lower spur threshold (ie more connectivity between blobs).

The idea that came up recently is: is it possible that the gap threshold is functionally getting defined by the difference between those two numbers? So rather than a gapThresh of 11, the 11/9 plot is applying a spatial kernel based on a gapThresh of [11-9 = 2], and a spurThresh of 11? I hope I've explained that decently well, please let me know if I can clarify! Thanks, Paige

[bi0m3trics]

@paigekouba I follow what you're saying, but at no point do I difference gap and spur... so that can't be it.

However, I did find where Brian Underwood provided the fragPatch.py code that was derived from the original patchMorph.py code, so I used it to strictly rework the whole thing. Have a look here In it you'll find a canopy height raster, and a script that is the strictest interpretation I could create for the original patchMorph.py code I had (which is in there too). I commented the shit out of the code, so it should be self-explanatory. Give it a run and let me know if it looks more like what you expected ( and what Danny and Jamie produced).

In make changes to explore helping you, I think I messed something up and now patchMorphSummary.R is acting weird so I just rewrote the hole thing. I'm also fairly certain it's faster, simpler, and more straight forward.

[paigekouba]

Hello again, I just wanted to let you know I have been waylaid by other projects but I will give this a try in the next week! I so appreciate the time you've dedicated to helping with my issue—I will update you with what I find out!

[bi0m3trics]

@paigekouba Any updates? I'd like to work on pushing the code I redid here if it seems to be working better., As I said, it's easier to read, simpler and actually faster.
Also, I'd like to get this issue closed if it's no longer being explored.

[paigekouba]

Hi Andrew, I have tried the new code and am still having the issue I mentioned, where the parameters don't seem to match the results. As you said the code is very clear and well commented, so I think I'll be able to get to the bottom of it. I'm in the finishing stretch of my dissertation and have had to focus on other projects, but I hope to return to this next month. I apologize for the slow progress!

[bi0m3trics]

I'm going to close this one for now (it can always be reopened), and I may migrate to the new code in version 1.0.3. In the end, I can't really be responsible for trying to match the results of someone else's code that's not available for comparison.