rCarto
commented
3 years ago Hello, st_within_distance() may be helpful here. https://r-spatial.github.io/sf/reference/geos_binary_pred.html
Closed k6adams closed 3 years ago
rCarto
commented
3 years ago Hello, st_within_distance() may be helpful here. https://r-spatial.github.io/sf/reference/geos_binary_pred.html
k6adams
commented
3 years ago Thanks @rCarto! That might be a good option. Still hoping I can’t figure out the bearing part of the equation. I am hoping to use some supervised models to estimate the radius of historical cyclone data. I suspect, but do not know, that the distance to land fall might be a meaningful predictor. Trouble is I don’t think the nearest piece of land matters as much as the nearest piece of land in a given direction. Like a cyclone might be closest to the Florida, but if is heading to the Carolinas and is still mostly over water, I am not sure I care that the closest piece of land is Florida.
edzer
commented
3 years ago Looks like this is not something you can do with sf.
k6adams
commented
3 years ago @edzer @rCarto I wanted to add my work-around incase anyone ever runs across this post and wants to do something similar.
library(tidyverse)
library(sf)
library(geosphere)
library(units)
library(lubridate)
## Create Lon Lat Vars
storm_tracks2 <- storm_tracks %>%
group_by(storm_id) %>%
arrange(date_time) %>%
mutate(
lat = map(geometry, 2) %>% unlist(),
lon = map(geometry, 1) %>% unlist(),
lat_lead = lead(lat),
lon_lead = lead(lon),
)
## Create Bearing Function
bearing_map <- function(geometry, lon_lead, lat_lead) {
bearing_accurate = bearing(
st_coordinates(geometry),
c(lon_lead, lat_lead)
) %>% units::set_units(degrees)
}
## Create Bearing to next point
storm_tracks3 <- storm_tracks2 %>%
group_by(storm_id) %>%
mutate(
bearing_accurate =
pmap(
.l = list(geometry=geometry, lon_lead=lon_lead, lat_lead=lat_lead),
.f = bearing_map
)
) %>%
ungroup()
## Create paths given a point and a specific bearing
gCB <- function(lon, lat, bearing_accurate) {
greatCircleBearing(
c(lon, lat),
bearing_accurate) %>%
st_linestring()
}
storm_tracks_gCB <- storm_tracks3 %>%
group_by(storm_id) %>%
mutate(
geometry=
pmap(
.l = list(lon=lon, lat=lat, bearing_accurate=bearing_accurate),
.f = gCB
) %>% st_as_sfc(crs=4326),
) %>%
st_as_sf(crs=4326) %>%
filter(!bearing_accurate== 0 %>% units::set_units(degrees))
## Could also just use the actual storm path
## I decided I prefer this
storm_paths <- storm_tracks3 %>%
group_by(storm_id) %>%
arrange(date_time) %>%
summarise(
elapsed_time = last(date_time) - first(date_time),
do_union = FALSE
) %>%
filter(elapsed_time > 0) %>%
st_cast('LINESTRING') %>%
mutate_if(is.numeric, list(~na_if(., Inf))) %>%
mutate_if(is.numeric, list(~na_if(., -Inf))) %>%
ungroup
## create map to measure point distance to storm
## north_america, is a shape with small land masses filtered out.
north_america_line <- north_america %>%
st_cast(to = 'POLYGON') %>%
st_cast(to = 'LINESTRING')
## find storm paths landfall point land.
landfall_points <- st_intersection(storm_paths, north_america_line)
st_crs(landfall_points)
ggplot() +
geom_sf(data = north_america, fill = NA) +
geom_sf(data = storm_paths) +
geom_sf(data = landfall_points, color=alpha("red", 0.2))
## Calculate distance to nearest landfall point
landfall_points2 <- landfall_points %>%
tibble %>%
rename(geometry_2 = geometry) %>%
mutate(
landfall = 1,
)
storm_tracks_landfall <- left_join(storm_tracks3, landfall_points2)
## Create Dist Function
dist_2_nearest <- function(geometry, geometry_2) {
st_nearest_points(geometry, geometry_2) %>% st_length
}
## Calculate Distance
storm_tracks_landfall_2 <- storm_tracks_landfall %>%
mutate(
dist_to_nearest_landfall_point =
pmap(
.l = list(geometry=geometry, geometry_2=geometry_2),
.f = dist_2_nearest
)
)
If you add the command that makes storm_tracks available, your example might be reproducible.
k6adams
commented
3 years ago If you add the command that makes
storm_tracksavailable, your example might be reproducible.
@edzer I will do that, but it may take me a week or two before I have time to post it. I vaguely know how to create a sample storm_tracks sf.
Should I include the north_america sf creation, or should I leave it as a placeholder?
edzer
commented
3 years ago That's up to you, but if you want to leave anything useful here, then leave a reprex.
BajczA475
commented
2 years ago I just found this thread--is there still no means by which to calculate distances between points and polygons in a specific bearing/direction?
k6adams
commented
2 years ago @edzer , @rCarto I know it has been a while, but here is my best attempt at a reprex. The data is being pulled from the web. I tried to use {datapasta} to embed the data in the code, but I had difficulty fixing the geometry column using that method.
@BajczA475 perhaps you might find this useful?
library(tidyverse)
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(geosphere)
library(units)
#> udunits database from C:/Users/Kevin.Adams3/OneDrive - USDA/Documents/R/win-library/4.1/units/share/udunits/udunits2.xml
library(lubridate)
#>
#> Attaching package: 'lubridate'
#> The following objects are masked from 'package:base':
#>
#> date, intersect, setdiff, union
library(curl)
#> Using libcurl 7.64.1 with Schannel
#>
#> Attaching package: 'curl'
#> The following object is masked from 'package:readr':
#>
#> parse_date
library(janitor)
#>
#> Attaching package: 'janitor'
#> The following objects are masked from 'package:stats':
#>
#> chisq.test, fisher.test
library(rnaturalearth)
library(rnaturalearthdata)
library(tigris)
#> To enable
#> caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.
library(nngeo)
## Download and prep data
ibtracs_shape_zip <- tempfile()
curl_download(
'https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs/v04r00/access/shapefile/IBTrACS.since1980.list.v04r00.points.zip',
ibtracs_shape_zip
)
ibtracs_shape_unzip <- tempfile()
unzip(
zipfile = ibtracs_shape_zip,
exdir = ibtracs_shape_unzip
)
remove(ibtracs_shape_zip)
ibtracs_shape <- list.files(
ibtracs_shape_unzip,
pattern = ".shp$",
full.names = TRUE
)
remove(ibtracs_shape_unzip)
ibtracs_4326 <- read_sf(ibtracs_shape) %>%
clean_names() %>%
select(
usa_atcfid,
usa_wind,
usa_pres,
name,
iso_time,
starts_with("usa_r64")
) %>%
filter(
usa_atcfid == "AL092020"
) %>%
mutate(
name = name %>% str_to_title(),
storm_name = str_c(
name,
usa_atcfid %>% str_trunc(4, side = "left", ellipsis = ""),
sep = ", "
) %>%
str_trim(),
iso_time = iso_time %>% as_datetime(),
) %>%
select(-name) %>%
mutate(
across(
starts_with("usa_r"),
as.numeric
)
) %>%
group_by(usa_atcfid) %>%
arrange(iso_time) %>%
mutate(
lat = map(geometry, 2) %>% unlist(),
lon = map(geometry, 1) %>% unlist(),
lat_lead = lead(lat),
lon_lead = lead(lon),
) %>%
ungroup()
## Create Bearing Function
bearing_map <- function(geometry, lon_lead, lat_lead) {
bearing_accurate = bearing(
st_coordinates(geometry),
c(lon_lead, lat_lead)
) %>% units::set_units(degrees)
}
## Create Bearing to next point
ibtracs_4326_2 <- ibtracs_4326 %>%
group_by(usa_atcfid) %>%
mutate(
bearing_accurate =
pmap(
.l = list(geometry=geometry, lon_lead=lon_lead, lat_lead=lat_lead),
.f = bearing_map
) %>% as.numeric()
) %>%
ungroup()
## Create paths given a point and a specific bearing
gCB <- function(lon, lat, bearing_accurate) {
greatCircleBearing(
c(lon, lat),
bearing_accurate,
n = 360
) %>%
st_linestring()
}
ibtracs_4326_gCB <- ibtracs_4326_2 %>%
drop_na(bearing_accurate) %>%
group_by(usa_atcfid) %>%
mutate(
geometry = case_when(
!is.na(lon_lead) & !is.na(lat_lead) ~
pmap(
.l = list(lon=lon, lat=lat, bearing_accurate=bearing_accurate),
.f = gCB
)
) %>% st_sfc(crs=4326)
) %>%
st_sf(crs = 4326) %>%
st_crop(xmin=-120, xmax=-40, ymin=10, ymax=60)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
ggplot()+
geom_sf(
data = ibtracs_4326_gCB %>% st_transform(5070)
) +
geom_sf(
data = ibtracs_4326_2 %>% st_transform(5070),
color ="blue",
size = 3
)
## Could also just use the actual storm path
## I decided I prefer this
ibtracs_4326_storm_paths <- ibtracs_4326_2 %>%
group_by(usa_atcfid) %>%
arrange(iso_time) %>%
summarise(
elapsed_time = last(iso_time) - first(iso_time),
do_union = FALSE
) %>%
filter(elapsed_time > 0) %>%
select(-elapsed_time) %>%
st_cast('LINESTRING') %>%
mutate_if(is.numeric, list(~na_if(., Inf))) %>%
mutate_if(is.numeric, list(~na_if(., -Inf))) %>%
ungroup
## create map to measure point distance to storm
## north_america, is a shape with small land masses filtered out.
usa_map_shape_5070 <- states(cb = TRUE) %>%
filter(
!(
STATEFP %in%
c(
"02",
"15",
"72",
"66",
"78",
"60",
"69",
"64",
"68",
"70",
"74",
"81",
"84",
"86",
"87",
"89",
"71",
"76",
"95",
"79")
)
) %>%
st_transform(5070) %>%
summarise(aland_sum = sum(ALAND)) %>%
select(geometry) %>%
st_cast("POLYGON") %>%
mutate(area = st_area(.)) %>%
arrange(desc(area)) %>%
slice(1) %>%
select(geometry)
#> | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
usa_map_shape_4326 <- st_transform(
usa_map_shape_5070,
crs = 4326
)
north_america_5070 <- ne_countries(
returnclass = "sf",
continent ="north america",
type = 'map_units',
scale=50
) %>%
st_transform(5070) %>%
mutate(
land_area = st_area(geometry)
)
north_america_5070_filter <- north_america_5070 %>%
mutate(across(where(is.character), str_trim)) %>%
filter(
!admin %in% c(
"Greenland",
#"Canada",
"United States of America")
) %>%
select(geometry) %>%
st_union()
ggplot() +
geom_sf(
data = north_america_5070_filter,
fill = NA
) +
geom_sf(
data = usa_map_shape_5070,
fill = NA
)
north_america_5070 <- st_union(usa_map_shape_5070, north_america_5070_filter) %>%
st_cast("POLYGON") %>%
mutate(area = st_area(.)) %>%
arrange(desc(area)) %>%
slice_head(n=27) %>%
st_union() %>%
st_remove_holes() %>%
st_simplify(dTolerance=1000, preserveTopology=T) %>%
st_as_sf()
north_america_4326 <- north_america_5070 %>%
st_transform(4326) %>%
st_as_sf()
ggplot() +
geom_sf(data = north_america_5070, fill = NA)
north_america_4326_line <- north_america_4326 %>%
st_cast(to = 'POLYGON') %>%
st_cast(to = 'LINESTRING')
## find storm paths landfall point land.
ibtracs_4326_landfall_points <- st_intersection(
ibtracs_4326_storm_paths,
north_america_4326_line
) %>%
group_by(
usa_atcfid
) %>%
summarise() %>%
ungroup()
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
st_crs(ibtracs_4326_landfall_points)
#> Coordinate Reference System:
#> User input: WGS 84
#> wkt:
#> GEOGCRS["WGS 84",
#> DATUM["World Geodetic System 1984",
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433]],
#> CS[ellipsoidal,2],
#> AXIS["latitude",north,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> AXIS["longitude",east,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> ID["EPSG",4326]]
ggplot() +
geom_sf(
data = north_america_4326,
fill = NA
) +
geom_sf(
data = ibtracs_4326_storm_paths
) +
geom_sf(
data = ibtracs_4326_landfall_points,
color = alpha("blue", 0.3),
size = 2
)
## Calculate distance to nearest landfall point
ibtracs_4326_landfall_points_2 <- ibtracs_4326_landfall_points %>%
rename(geometry_2 = geometry)
ibtracs_4326_landfall <- left_join(
ibtracs_4326_2,
ibtracs_4326_landfall_points_2 %>%
tibble() %>%
select(
geometry_2,
usa_atcfid
)
) %>%
st_sf(crs=4326)
#> Joining, by = "usa_atcfid"
## Create Dist Function
dist_2_nearest <- function(geometry, geometry_2) {
st_nearest_points(geometry, geometry_2) %>% st_length
}
## Calculate Distance
ibtracs_4326_landfall_2 <- ibtracs_4326_landfall %>%
mutate(
dist_to_nearest_landfall_point =
pmap(
.l = list(geometry=geometry, geometry_2=geometry_2),
.f = dist_2_nearest
),
) Created on 2022-02-24 by the reprex package (v2.0.1)
BajczA475
commented
2 years ago TY for thinking of me. I'll be excited to take a look at this, although for the moment this post details an approach that works for me: https://stackoverflow.com/questions/71240972/how-can-i-draw-a-line-from-a-point-to-a-polygon-edge-and-then-get-the-lines-len/71241675?noredirect=1#comment125939777_71241675
I am new to R, GitHub, and GIS in general, so I apologize a head of time if this is the wrong place to ask this question.
I am wondering if there is a way calculate the distance between a point and a polygon with a specified bearing?
Specifically, I am looking to calculate a tropical cyclone's distance to land given a specific bearing.
I know how to calculate a cyclone's distance to land using st_distance(). I know how to calculate a cyclone's bearing using two points and the bearing() function from the
geospherepackage. (thank you to @edzer).Does anyone have an idea how I might mash these functions together? Also, my hope would be set a max distance, so in the case where the specified bearing did not lead to polygon, the function would return an NA value.