spatstat.geom::owin() dependency
Closed this issue · 5 comments
In makeLines, dependencies on spatstat.geom tend to cause errors. There have been several cases where a cutout polygon is not seen by spatstat.geom::owin(). The latest is the meadowdale cutout kml.
spatstat.geom::owin() fails to respect the cutout:
However, this can be fixed many different ways but no obvious reason as to why solutions work. One could add a buffer:
library(dplyr)
library(pronghornLT)
library(sf)
# data
sPoly <- huntAreas %>%
filter(HERDNAME == "Meadowdale")
xPoly <- st_read("C:/Users/garcatlin/Desktop/Meadowdale.LT_Exlusion.Polygon.kml")
# buffer
xPoly <- st_transform(xPoly, st_crs(sPoly))
xPoly <- st_buffer(xPoly, 100)
One could take the difference between the huntarea and the exclusion and union them back together:
# data
sPoly <- huntAreas %>%
filter(HERDNAME == "Meadowdale")
xPoly <- st_read("C:/Users/garcatlin/Desktop/Meadowdale.LT_Exlusion.Polygon.kml")
# difference
diff <- st_difference(xPoly, sPoly)
xPoly <- st_union(diff, xPoly)
Or finally, disabling spatstat's attempt at fixing polygons also works (though this spits out plenty of warning messages in the console):
library(spatstat.geom)
spatstat.options(fixpolygons=FALSE)
At any rate, a revamp of makeLines and associated helper functions that would bypass reliance on spatstat.geom might be in order
@jcarlis3 Will you check with my open branch? I think this solves Lee's issue and also allows for our not passing an xpoly to makeLines.
Here's what I was testing with:
library(pronghornLT)
library(dplyr)
library(sf)
library(mapview)
# herd units
s = herdUnits %>%
filter(HERDNAME %in% c("Badwater", "North Natrona"))
# make lines
l = makeLines(s, totalLengthKm = 2800)
mapview(s) + mapview(l$lines, color = "red")
# cutouts
x = st_read("~/Documents/Github/pronghornLT_app/Data/GIS/Badwater/BW Ant UnOc Habitat 2016 UTM13.shp")
x2 = st_read("~/Documents/Github/pronghornLT_app/Data/GIS/North Natrona/NorthNatrona_HU_ToRemove.shp")
x3 = st_union(st_geometry(x), st_geometry(x2))
# map
mapview(s) +
mapview(x3, col.regions = "red")
# use cutPoly on x3
diff = cutPoly(s, x3)
# map
mapview(diff)
# pass to calc line length
line_rec = calcLineLength(
diff,
15000
)
# make lines
l = makeLines(diff,
totalLengthKm = line_rec)
# map
mapview(diff) +
mapview(l$lines, color = "red")
# use cutPoly on x only
diff = cutPoly(s, x)
# map
mapview(diff)
# pass to calc line length
line_rec = calcLineLength(
diff,
15000
)
# make lines
l = makeLines(diff,
totalLengthKm = line_rec,
considerHuntArea = T)
# map
mapview(diff) +
mapview(l$lines, color = "red")
Here's Lee's from this afternoon:
library(pronghornLT)
library(dplyr)
library(sf)
library(mapview)
# herd units
s = herdUnits %>%
filter(HERDNAME %in% c("Medicine Bow"))
# cutouts
x = st_read("~/Desktop/out polygons.kml")
# map
mapview(s) +
mapview(x, col.regions = "red")
# use cutPoly on x3
diff = cutPoly(s, x)
# map
mapview(diff)
# pass to calc line length
line_rec = calcLineLength(
diff,
15000
)
# make lines
l = makeLines(diff,
totalLengthKm = line_rec)
# plot
mapview(s) +
mapview(x, col.regions = "red") +
mapview(l$lines)
@grcatlin , nicely done! I tested on Sublette on the original pathological cutout shp. The only issue I see is that makeLines appears to no longer be able to break the lines at the hunt area boundaries.
# Test new cutPolys with Sublette
library(pronghornLT)
library(dplyr)
library(sf)
library(mapview)
# herd units
s = herdUnits %>%
filter(HERDNAME %in% c("Sublette"))
# cutouts
x <- st_read("C:/Users/jadcarlisle/My Drive/Pronghorn/PronghornLT/TransectLayouts/Sublette/problematic/2023 Sublette LT Cutout Polygon.shp")
# map
mapview(s) +
mapview(x, col.regions = "red")
# use cutPoly on x3
diff = cutPoly(s, x)
# map
mapview(diff)
# pass to calc line length
line_rec = calcLineLength(
diff,
N = 20e3
)
# make lines
l = makeLines(diff,
totalLengthKm = line_rec,
considerHuntArea = TRUE)
ha <- huntAreas %>%
filter(HERDNAME == s$HERDNAME)
# plot
mapview(s) +
mapview(ha, alpha.regions = 0) +
mapview(x, col.regions = "red") +
mapview(l$lines)
@jcarlis3 What happened here is that previously I had relied on makeLines being based an uncut polygon. I used some sf maths to get the outer borders of that polygon and do slicing on huntArea lines to calculate where the transect lines intersected inner boundaries. Since we are no longer passing a whole polygon, this got a lot more complicated.
What I tried yesterday was taking the centroids of the differenced polygon (which requires it to not be unioned) to see which hunt areas the centroid intersect so I could back-calculate the whole herd unit, but this didn't work so well as you saw. Unless there's a better way you can think of, I've basically made it so users will be required to use the huntAreas data within the package. This allows us to do some filtering so we can extract out the names of the huntAreas in sPoly and use those to create the whole polygon to calculate hunt area intersects.
TL;DR: sPoly must be originally filtered from pronghornLT::huntAreas to allow for us to calculate intersections. Not sure if this is a bad thing, but seeing as how the app no longer takes uploads for herdunit/huntarea shapefiles anyway, I think it makes a little sense. Do let me know if there's a better way you can think of though!
I also added some mapping functions as that will allow us to reduce code in the app.
library(pronghornLT)
library(dplyr)
library(sf)
library(mapview)
# badwater & nn -----------------------------------------------------------
# herd units
s = huntAreas %>%
filter(HERDNAME %in% c("Badwater", "North Natrona"))
# cutouts
x = st_read("~/Documents/Github/pronghornLT_app/Data/GIS/Badwater/BW Ant UnOc Habitat 2016 UTM13.shp")
x2 = st_read("~/Documents/Github/pronghornLT_app/Data/GIS/North Natrona/NorthNatrona_HU_ToRemove.shp")
x3 = st_union(st_geometry(x), st_geometry(x2))
# use cutPoly on x3
diff = cutPoly(s, x3)
# pass to calc line length
line_rec = calcLineLength(
diff,
15000
)
# make lines
l = makeLines(diff,
totalLengthKm = line_rec)
# map
mapHerdWithLines(s, x3, l$lines)
# med bow -----------------------------------------------------------------
# herd units
s = huntAreas %>%
filter(HERDNAME %in% c("Medicine Bow"))
# cutouts
x = st_read("~/Desktop/out polygons.kml")
# use cutPoly on x3
diff = cutPoly(s, x)
# pass to calc line length
line_rec = calcLineLength(
diff,
15000
)
# make lines
l = makeLines(diff,
totalLengthKm = line_rec)
# map
mapHerdWithLines(s, x, l$lines)
# sublette ----------------------------------------------------------------
# herd units
s = huntAreas %>%
filter(HERDNAME %in% c("Sublette"))
# cutouts
x <- st_read("~/Documents/Github/pronghornLT_app/Data/GIS/Sublette/2023 Sublette LT Cutout Polygon.shp")
# use cutPoly on x3
diff = cutPoly(s, x)
# pass to calc line length
line_rec = calcLineLength(
diff,
15000
)
# make lines
l = makeLines(diff,
totalLengthKm = line_rec)
# map
mapHerdWithLines(s, x, l$lines)
# other maps
mapHerd(s)
mapHerdWithCutouts(s, x)
# herd without cutouts and lines ------------------------------------------
s = huntAreas %>%
filter(HERDNAME == "Leiter")
# pass to calc line length
line_rec = calcLineLength(
s,
15000
)
# make lines
l = makeLines(s,
totalLengthKm = line_rec)
# map
mapHerdWithLines(herdPoly = s, lines = l$lines)
@jcarlis3 Also I just updated the map functions to be one instead of multiple so it'll look like this now:
# herd units
s = huntAreas %>%
filter(HERDNAME %in% c("Badwater", "North Natrona"))
# cutouts
x = st_read("~/Documents/Github/pronghornLT_app/Data/GIS/Badwater/BW Ant UnOc Habitat 2016 UTM13.shp")
x2 = st_read("~/Documents/Github/pronghornLT_app/Data/GIS/North Natrona/NorthNatrona_HU_ToRemove.shp")
x3 = st_union(st_geometry(x), st_geometry(x2))
x3 = st_as_sf(x3)
# use cutPoly on x3
diff = cutPoly(s, x3)
# pass to calc line length
line_rec = calcLineLength(
diff,
15000
)
# make lines
l = makeLines(diff,
totalLengthKm = line_rec)
# map
mapLT(s, x, l$lines)
mapLT(s)
mapLT(s, x)