BenioffOceanInitiative / bbnj-scripts

migrating bbnj:inst/scripts to save space in R package bbnj
https://BenioffOceanInitiative.github.io/bbnj-scripts
MIT License
0 stars 0 forks source link

error when updating scenario_overlays #1

Closed mvisalli closed 5 years ago

mvisalli commented 5 years ago

Hi @bbest,

I'm attempting to rerun scenario_overlays.Rmd with the new prioritir s4 solution (including Antarctica) in order to get updated % overlap of s4 with each IHO sea for manuscript.

When I get to the for loop on line 89:

for (layer in lyrs_clip$layer){ # layer = lyrs_clip$layer[1]
    lyr <- layers %>% filter(layer==!!layer)
    read_sf(lyr$raw_shp) %>%
      st_intersection(p_abnj) %>% 
      mutate(
        area_km2 = st_area(geometry) %>%
          units::set_units(km^2)) %>% 
      filter(
        st_geometry_type(geometry) != "GEOMETRYCOLLECTION") %>% # to remove points
      write_sf(lyr$hs_shp)
  }

I receive the following error:

although coordinates are longitude/latitude, st_intersection assumes that they are planar

Error in CPL_geos_op2(op, st_geometry(x), st_geometry(y)) : 
  Evaluation error: TopologyException: Input geom 0 is invalid: Ring Self-intersection at or near point 3.5156249243664486 43.274444391192986 at 3.5156249243664486 43.274444391192986.

I think it may be having trouble with p_abnj on line 93 -- perhaps since this .Rmd file is now in the repo bbnj-scripts instead of bbnj? But I'm not sure how to resolve this. Help??

mvisalli commented 5 years ago

On a related note, would it be possible to use this same script to evaluate s4 boundaries with EEZs to get # of EEZs that touch the s4 solution? Perhaps it could be similar to how it evaluates overlap with WDPA MPAs (that are not in hs) such as:

Screen Shot 2019-10-30 at 11 11 06 AM

Do you have an updated eez shapefile with Antarctica removed? I can't seem to find it in the bbnj repo

mvisalli commented 5 years ago

@bbest do we need to be concerned about using st_intersection on a sphere? see this github issue

sf is informed that these are long/lat Earth coordinates, and the arctic polygon is supposed to include the pole, but st_intersects says it doesn't, which is wrong. It's wrong because it assumes coordinates are planar, and you are warned about that assumption.

(reality is somewhat more complicated: every polygon on the sphere divides the sphere into two parts, depending on the ring direction; this would have to be taken into account)

I hope to extend sf with routines that do proper intersections on the sphere; R package s2 from Ege @rubak seems to be a good candidate for this.

See also: r-spatial.org spherical geometry

bbest commented 5 years ago

Hi @mvisalli,

I think I was able successfully convert the GEOMETRYCOLLECTION of North Pacific in ihor.shp into MULTIPOLYGON by using st_collection_extract("POLYGON") then group_by(...) %>% summarize(). Need to jump on another call but wanted to share good news and will check in fixed outputs this afternoon.

Ben

mvisalli commented 5 years ago

Thanks @bbest I'll keep an eye out for updated results & will plug into manuscript when ready

bbest commented 5 years ago

Hi @mvisalli,

All set! Check out North Pacific Ocean now in the IHO seas revised (ihor):

https://ecoquants.com/bbnj-scripts/scenario_overlays.html#ihor

image