Fetch in a bay

[lzhu5]

Hi, I am trying to get a fetch in Delaware Bay in the U.S. I followed one test case (Philippine case) from the manual. I got the coastline using the map package as follows: usa.df = ggplot2::map_data("usa", region = "main")

My error message is: Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false, : TopologyException: Input geom 1 is invalid: Self-intersection at or near point -3709170.8294415493 5439551.5041889455 at -3709170.8294415493 5439551.5041889455

Could you please enlighten me how to fix this problem? Thanks.

[blasee]

Hi @lzhu5,

I am unable to reproduce this error, and this error does not come from the fetchR package. Can you please show me the output of sessionInfo()? Are all of your packages and R version up to date?

I am able to calculate fetch at a random site within Delaware Bay using the following R code (UPDATED):

library(fetchR)
library(ggplot2)
library(sp)

delaware_bay = map_data("state", region = c("delaware", "new jersey"))

ggplot(delaware_bay, aes(long, lat, group = group)) +
  geom_point(aes(-75.3, 39.2), shape = 3, colour = "red", size = 3) +
  geom_polygon() +
  coord_quickmap()

# Uses the NAD83 (2011) / Delaware EPSG 6435 map projection.
# 
# Create a data frame of the locations at which to calculate fetch
fetch_locs = data.frame(long = -75.3, lat = 39.2)

# Convert to a spatial points object.
# Note the CRS of the locations are lat/lon (epsg 4324). We are then 
# transforming the CRS to EPSG 6435!
delaware_loc = spTransform(SpatialPoints(fetch_locs, CRS("+init=epsg:4324")),
                           CRS("+init=epsg:6435"))

# Create a spatial polygon object for the coastline
delaware_poly = Polygon(subset(delaware_bay, region == "delaware", 1:2))
new_jersey_poly = Polygon(subset(delaware_bay, region == "new jersey", 1:2))

# Note that the CRS is also in lat/lon
delaware_bay_sp = SpatialPolygons(list(Polygons(list(delaware_poly, new_jersey_poly), 
                                                ID = "delaware bay")),
                                  proj4string = CRS("+init=epsg:4324"))

# Make sure we have the correct locations and polygons in lat/lon space
plot(delaware_bay_sp, col = "lightgrey", border = NA)
points(fetch_locs, col = 2, pch = 3, cex = 2)
axis(1)
axis(2)

# Make sure we have the correct locations in the projected space
plot(spTransform(delaware_bay_sp, CRS("+init=epsg:6435")),
     col = "lightgrey", border = NA)
points(delaware_loc, col = 2, pch = 3, cex = 2)
axis(1)
axis(2)

# Calculate fetch
delaware_bay_fetch = fetch(delaware_bay_sp, delaware_loc)

delaware_bay_fetch

Is projected    : TRUE
Max distance    : 300 km
Directions  : 36
Sites       : 1

       North East South West Average
Site 1  17.3   23  88.3   10    34.6

plot(delaware_bay_fetch, add = TRUE)

# Get a dataframe of the distances for each vector (in km)
as(delaware_bay_fetch, "data.frame")

Session Information

sessionInfo()
R version 3.5.2 (2018-12-20)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)

Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.7.0
LAPACK: /usr/lib/lapack/liblapack.so.3.7.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] maps_3.3.0    sp_1.3-1      ggplot2_3.1.0 fetchR_2.1-1 

loaded via a namespace (and not attached):
 [1] splines_3.5.2       Formula_1.2-3       assertthat_0.2.0    latticeExtra_0.6-28 pixmap_0.4-11      
 [6] pillar_1.3.1        backports_1.1.3     lattice_0.20-38     glue_1.3.0          digest_0.6.18      
[11] RColorBrewer_1.1-2  checkmate_1.8.5     colorspace_1.3-2    htmltools_0.3.6     Matrix_1.2-15      
[16] plyr_1.8.4          XML_3.98-1.16       aqp_1.17            pkgconfig_2.0.2     raster_2.8-4       
[21] purrr_0.2.5         scales_1.0.0        intervals_0.15.1    gstat_1.1-6         htmlTable_1.13.1   
[26] tibble_2.0.0        RSAGA_1.3.0         withr_2.1.2         nnet_7.3-12         lazyeval_0.2.1     
[31] survival_2.42-4     magrittr_1.5        crayon_1.3.4        MASS_7.3-51.1       xts_0.11-2         
[36] foreign_0.8-71      class_7.3-14        FNN_1.1.2.2         tools_3.5.2         dismo_1.1-4        
[41] shapefiles_0.7      data.table_1.11.8   stringr_1.3.1       munsell_0.5.0       cluster_2.0.7-1    
[46] plotrix_3.7-4       bindrcpp_0.2.2      colorRamps_2.3      compiler_3.5.2      e1071_1.7-0        
[51] spacetime_1.2-2     rlang_0.3.1         classInt_0.3-1      grid_3.5.2          rstudioapi_0.8     
[56] htmlwidgets_1.3     base64enc_0.1-3     labeling_0.3        gtable_0.2.0        codetools_0.2-15   
[61] reshape_0.8.8       R6_2.3.0            gridExtra_2.3       zoo_1.8-4           knitr_1.21         
[66] dplyr_0.7.8         rgdal_1.3-6         rgeos_0.4-2         bindr_0.1.1         Hmisc_4.1-1        
[71] stringi_1.2.4       Rcpp_1.0.0          plotKML_0.5-8       rpart_4.1-13        acepack_1.4.1      
[76] tidyselect_0.2.5    xfun_0.4

[lzhu5]

Hi Blake,

Thanks a lot for your reply. My sessionInf() outputs are very much like yours.

After comparing your code and mine, I found that the map_data caused the problem. If I switch to map_data("state", region = c("delaware", "new jersey")) as you suggested, I can calculate the fetch.

Thanks again for your help.

[lzhu5]

Hi Blake,

I revisited the fetch results that you posted. Why does the north and west fetch = 0? From the vector plot, their lengths are not zero.

[blasee]

Hi @lzhu5,

Sorry for the late reply.

You are correct, there shouldn't have been a fetch of 0, and it turns out I had a mistake in the code I gave you (sorry...). I've updated the code in my above comment. The mistake was that I had used lat/lon coordinates to define a point in the projected space. You will see that the code now correctly defines the point in lat/lon space, and then spTranform()'s into the projected space.

I hope this helps.

Cheers, Blake

[lzhu5]

Hi Blake,

Thanks a lot for your help!

Best, Ling