Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

st_intersection() gives Error in CPL_nary_intersection(x) : GEOS exception #1668

Closed
MHaringa opened this issue May 3, 2021 · 18 comments
Closed

Comments

@MHaringa
Copy link

MHaringa commented May 3, 2021

For example, I have the following data set:

df <- structure(list(lat = c(53.218461547511, 53.2354543814792, 53.2219368310548,
53.2142084420887, 53.2051969857393), lon = c(6.57022923569802,
6.55005158685787, 6.57015231400986, 6.5633788874466, 6.5692156138519
)), row.names = c(NA, -5L), class = c("tbl_df", "tbl", "data.frame"
))

Then I create a circle around each point:

circle_df <- df %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
st_transform(3035) %>%
st_buffer(dist = units::set_units(50, "kilometers"))

Applying st_intersection() gives an error:

circle_df %>%
st_intersection()

Error in CPL_nary_intersection(x) : GEOS exception

If I only use the first four rows no error message is returned:

circle_df[1:4,] %>%
st_intersection()

I don't get why an error message is returned here.

@edzer
Copy link
Member

edzer commented May 3, 2021

Because GEOS runs into an exception. You created a very difficult problem for it, with data like this:
x

@edzer
Copy link
Member

edzer commented May 3, 2021

This reminds me btw that #1535 still needs to be done, and might be needed to get the new geom engine to do its work.

edzer added a commit that referenced this issue May 3, 2021
@edzer
Copy link
Member

edzer commented May 3, 2021

You may want to try sth along the lines of

circle_df %>% st_set_precision(1) %>% st_intersection()

with this patch, that seems to work better than without.

@dbaston
Copy link
Contributor

dbaston commented May 3, 2021

Does this fail with the new overlay engine (GEOS >= 3.9) ?

@edzer
Copy link
Member

edzer commented May 3, 2021

Using branch geos_precision, without setting precision I see errors from GEOS 3.9.0:

library(sf)
# Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(dplyr)

# Attaching package: ‘dplyr’

# The following objects are masked from ‘package:stats’:

#     filter, lag

# The following objects are masked from ‘package:base’:

#     intersect, setdiff, setequal, union

df <- structure(list(lat = c(53.218461547511, 53.2354543814792, 53.2219368310548,
53.2142084420887, 53.2051969857393), lon = c(6.57022923569802,
6.55005158685787, 6.57015231400986, 6.5633788874466, 6.5692156138519
)), row.names = c(NA, -5L), class = c("tbl_df", "tbl", "data.frame"
))

#Then I create a circle around each point:

circle_df <- df %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
st_transform(3035) %>%
st_buffer(dist = units::set_units(50, "kilometers"))

#Applying st_intersection() gives an error:

#y <- circle_df %>% st_set_precision(1) %>% st_intersection()
circle_df %>% st_intersection()
# CBR: result (after common-bits addition) is INVALID: Self-intersection at or near point 4042052.512894101 3353674.4936534893 (4042052.5128941009752 3353674.4936534892768)
# <A>
# POLYGON ((4050058.9913444779813290 3324151.4958460354246199, 4048691.2495525269769132 3326383.4475967870093882, 4047442.1935323304496706 3328683.9226098093204200, 4046315.2468596189282835 3331046.6154429968446493, 4045313.4984168889932334 3333465.0501195220276713, 4044439.6939269914291799 3335932.5978780393488705, 4043696.2284272955730557 3338442.4953416609205306, 4043085.1397050586529076 3340987.8630558988079429, 4042608.1027119918726385 3343561.7243447750806808, 4042266.4249733351171017 3346157.0244334042072296, 4042061.0430040201172233 3348766.6497846394777298, 4041992.5197417489252985 3351383.4475967870093882, 4042052.5128941009752452 3353674.4936534911394119, 4042052.5128941009752452 3353674.4936534892767668, 4042052.5128941009752452 3353674.4936534948647022, 4042253.1055478826165199 3356223.2649810146540403, 4042594.7832865393720567 3358818.5650696437805891, 4043071.8202796061523259 3361392.4263585200533271, 4043682.9090018430724740 3363937.7940727584064007, 4044426.3745015384629369 3366447.6915363795123994, 4045300.1789914364926517 3368915.2392948972992599, 4046301.9274341664277017 3371333.6739714220166206, 4047428.8741068779490888 3373696.3668046095408499, 4048677.9301270744763315 3375996.8418176323175430, 4050045.6719190250150859 3378228.7935683834366500, 4051528.3505975487641990 3380386.1044322559610009, 4053121.9022434479556978 3382462.8613701239228249, 4054821.9590424266643822 3384453.3721355749294162, 4056623.8612569686956704 3386352.1808769595809281, 4058084.9531613294966519 3387738.7053229119628668, 4060033.5240568132139742 3389402.9420890007168055, 4062110.2809946816414595 3390996.4937348994426429, 4064267.5918585537001491 3392479.1724134231917560, 4066499.5436093052849174 3393846.9142053741961718, 4068800.0186223275959492 3395095.9702255707234144, 4071162.7114555151201785 3396222.9168982822448015, 4073581.1461320403032005 3397224.6653410121798515, 4076048.6938905576243997 3398098.4698309097439051, 4078558.5913541791960597 3398841.9353306056000292, 4081103.9590684170834720 3399453.0240528425201774, 4083677.8203572933562100 3399930.0610459093004465, 4086273.1204459224827588 3400271.7387845660559833, 4088882.7457971577532589 3400477.1207538810558617, 4091499.5436093052849174 3400545.6440161522477865, 4094116.3414214523509145 3400477.1207538810558617, 4096725.9667726876214147 3400271.7387845660559833, 4099321.2668613167479634 3399930.0610459093004465, 4101895.1281501930207014 3399453.0240528425201774, 4104440.4958644309081137 3398841.9353306056000292, 4106950.3933280524797738 3398098.4698309102095664, 4109417.9410865702666342 3397224.6653410126455128, 4111836.3757630949839950 3396222.9168982822448015, 4114199.0685962825082242 3395095.9702255707234144, 4116499.5436093052849174 3393846.9142053741961718, 4118731.4953600564040244 3392479.1724134236574173, 4120888.8062239289283752 3390996.4937348999083042, 4122965.5631617968901992 3389402.9420890011824667, 4124956.0739272478967905 3387702.8852900220081210, 4126854.8826686325483024 3385900.9830754799768329, 4128656.7848831750452518 3384002.1743340953253210, 4130356.8416821537539363 3382011.6635686443187296, 4131950.3933280524797738 3379934.9066307763569057, 4133433.0720065762288868 3377777.5957669038325548, 4134800.8137985272333026 3375545.6440161527134478, 4136049.8698187237605453 3373245.1690031299367547, 4137176.8164914352819324 3370882.4761699424125254, 4138178.5649341652169824 3368464.0414934176951647, 4138663.3753943806514144 3367094.9789877417497337, 4139020.0694404873065650 3365890.8000689102336764, 4139631.1581627242267132 3363345.4323546718806028, 4140108.1951557910069823 3360771.5710657956078649, 4140449.8728944477625191 3358176.2709771664813161, 4140655.2548637627623975 3355566.6456259312108159, 4140723.7781260339543223 3352949.8478137836791575, 4140655.2548637627623975 3350333.0500016366131604, 4140449.8728944477625191 3347723.4246504008769989, 4140108.1951557910069823 3345128.1245617722161114, 4139631.1581627242267132 3342554.2632728959433734, 4139020.0694404873065650 3340008.8955586575902998, 4138276.6039407914504409 3337498.9980950364843011, 4137402.7994508938863873 3335031.4503365186974406, 4136401.0510081639513373 3332613.0156599935144186, 4135274.1043354524299502 3330250.3228268064558506, 4134025.0483152559027076 3327949.8478137836791575, 4132657.3065233053639531 3325717.8960630325600505, 4131174.6278447811491787 3323560.5851991600356996, 4129581.0761988824233413 3321483.8282612916082144, 4127881.0193999037146568 3319493.3174958406016231, 4126079.1171853612177074 3317594.5087544564157724, 4124180.3084439770318568 3315792.6065399139188230, 4122189.7976785260252655 3314092.5497409352101386, 4120113.0407406575977802 3312498.9980950364843011, 4117955.7298767855390906 3311016.3194165127351880, 4115723.7781260339543223 3309648.5776245617307723, 4113423.3031130111776292 3308399.5216043652035296, 4111060.6102798241190612 3307272.5749316536821425, 4108642.1756032989360392 3306270.8264889237470925, 4106174.6278447816148400 3305397.0219990261830389, 4103664.7303811600431800 3304653.5564993303269148, 4101119.3626669221557677 3304042.4677770934067667, 4098545.5013780454173684 3303565.4307840266264975, 4095950.2012894167564809 3303223.7530453698709607, 4093340.5759381810203195 3303018.3710760548710823, 4090723.7781260339543223 3302949.8478137836791575, 4088106.9803138868883252 3303018.3710760548710823, 4085497.3549626511521637 3303223.7530453698709607, 4082902.0548740224912763 3303565.4307840266264975, 4080328.1935851462185383 3304042.4677770934067667, 4077782.8258709078654647 3304653.5564993303269148, 4075272.9284072867594659 3305397.0219990261830389, 4072805.3806487689726055 3306270.8264889237470925, 4070386.9459722442552447 3307272.5749316536821425, 4068024.2531390567310154 3308399.5216043652035296, 4065723.7781260339543223 3309648.5776245617307723, 4063491.8263752828352153 3311016.3194165122695267, 4061334.5155114103108644 3312498.9980950360186398, 4059257.7585735423490405 3314092.5497409352101386, 4057267.2478080913424492 3315792.6065399139188230, 4055368.4390667066909373 3317594.5087544564157724, 4053566.5368521641939878 3319493.3174958406016231, 4052524.1232459573075175 3320713.8270143778063357, 4051541.6700230017304420 3321994.1849821633659303, 4050058.9913444779813290 3324151.4958460354246199))
# </A>
# CBR: result (after common-bits addition) is INVALID: Self-intersection at or near point 4048532.4282185733 3376090.9350878298 (4048532.4282185733318 3376090.9350878298283)
# <A>
# POLYGON ((4042052.5128941014409065 3353674.4936534985899925, 4042061.0430040201172233 3354000.2454089340753853, 4042266.4249733351171017 3356609.8707601693458855, 4042608.1027119918726385 3359205.1708487984724343, 4043085.1397050586529076 3361779.0321376747451723, 4043696.2284272955730557 3364324.3998519130982459, 4044439.6939269909635186 3366834.2973155342042446, 4045313.4984168889932334 3369301.8450740519911051, 4046315.2468596189282835 3371720.2797505767084658, 4047442.1935323304496706 3374082.9725837642326951, 4048532.4282185742631555 3376090.9350878316909075, 4048532.4282185733318329 3376090.9350878298282623, 4048532.4282185742631555 3376090.9350878316909075, 4049566.0152120338752866 3377777.5957669033668935, 4051048.6938905576243997 3379934.9066307758912444, 4052642.2455364568158984 3382011.6635686438530684, 4054342.3023354355245829 3384002.1743340948596597, 4056144.2045499775558710 3385900.9830754795111716, 4058043.0132913622073829 3387702.8852900220081210, 4058084.9531613290309906 3387738.7053229114972055, 4056623.8612569686956704 3386352.1808769595809281, 4054821.9590424266643822 3384453.3721355749294162, 4053121.9022434479556978 3382462.8613701239228249, 4051528.3505975487641990 3380386.1044322559610009, 4050045.6719190250150859 3378228.7935683834366500, 4048677.9301270744763315 3375996.8418176323175430, 4047428.8741068779490888 3373696.3668046095408499, 4046301.9274341664277017 3371333.6739714220166206, 4045300.1789914364926517 3368915.2392948972992599, 4044426.3745015384629369 3366447.6915363795123994, 4043682.9090018430724740 3363937.7940727584064007, 4043071.8202796061523259 3361392.4263585200533271, 4042594.7832865393720567 3358818.5650696437805891, 4042253.1055478826165199 3356223.2649810146540403, 4042052.5128941014409065 3353674.4936534985899925))
# </A>
# CBR: result (after common-bits addition) is INVALID: Self-intersection at or near point 4048532.4282185747 3376090.9350878326 (4048532.4282185747288 3376090.9350878326222)
# <A>
# GEOMETRYCOLLECTION (LINESTRING (4134959.6351324799470603 3325838.1565251075662673, 4134800.8137985272333026 3325545.6440161522477865), POLYGON ((4042052.5128941009752452 3353674.4936534911394119, 4042047.7235785676166415 3353613.6396297793835402, 4042253.1055478826165199 3356223.2649810146540403, 4042052.5128941009752452 3353674.4936534911394119)), POLYGON ((4049331.1645958102308214 3324930.8600477352738380, 4050045.6719190254807472 3323764.8900668807327747, 4051528.3505975492298603 3321607.5792030086740851, 4053121.9022434479556978 3319530.8222651402465999, 4054821.9590424266643822 3317540.3114996892400086, 4056623.8612569691613317 3315641.5027583050541580, 4058522.6699983538128436 3313839.6005437625572085, 4060513.1807638048194349 3312139.5437447838485241, 4062589.9377016727812588 3310545.9920988846570253, 4064747.2485655453056097 3309063.3134203609079123, 4066979.2003162964247167 3307695.5716284103691578, 4069279.6753293192014098 3306446.5156082138419151, 4071642.3681625067256391 3305319.5689355023205280, 4074060.8028390314429998 3304317.8204927723854780, 4076528.3505975492298603 3303444.0160028748214245, 4079038.2480611703358591 3302700.5505031789653003, 4081583.6157754086889327 3302089.4617809420451522, 4084157.4770642849616706 3301612.4247878752648830, 4086752.7771529136225581 3301270.7470492185093462, 4089362.4025041493587196 3301065.3650799035094678, 4091979.2003162964247167 3300996.8418176323175430, 4094595.9981284434907138 3301065.3650799035094678, 4097205.6234796792268753 3301270.7470492185093462, 4099800.9235683078877628 3301612.4247878752648830, 4102374.7848571846261621 3302089.4617809420451522, 4104920.1525714225135744 3302700.5505031789653003, 4107430.0500350440852344 3303444.0160028748214245, 4109897.5977935614064336 3304317.8204927723854780, 4112316.0324700865894556 3305319.5689355023205280, 4114678.7253032736480236 3306446.5156082138419151, 4116979.2003162964247167 3307695.5716284103691578, 4119211.1520670480094850 3309063.3134203609079123, 4121368.4629309200681746 3310545.9920988851226866, 4123445.2198687884956598 3312139.5437447838485241, 4125393.7907642694190145 3313803.7805108702741563, 4124956.0739272483624518 3313388.4027422824874520, 4122965.5631617973558605 3311688.3459433037787676, 4120888.8062239289283752 3310094.7942974050529301, 4118731.4953600568696856 3308612.1156188808381557, 4116499.5436093052849174 3307244.3738269302994013, 4114199.0685962825082242 3305995.3178067337721586, 4111836.3757630954496562 3304868.3711340222507715, 4109417.9410865702666342 3303866.6226912923157215, 4106950.3933280529454350 3302992.8182013947516680, 4104440.4958644313737750 3302249.3527016988955438, 4101895.1281501934863627 3301638.2639794619753957, 4099321.2668613167479634 3301161.2269863951951265, 4096725.9667726880870759 3300819.5492477384395897, 4094116.3414214523509145 3300614.1672784234397113, 4091499.5436093052849174 3300545.6440161522477865, 4088882.7457971582189202 3300614.1672784234397113, 4086273.1204459224827588 3300819.5492477384395897, 4083677.8203572938218713 3301161.2269863951951265, 4081103.9590684175491333 3301638.2639794619753957, 4078558.5913541791960597 3302249.3527016988955438, 4076048.6938905580900609 3302992.8182013947516680, 4073581.1461320403032005 3303866.6226912923157215, 4071162.7114555155858397 3304868.3711340222507715, 4068800.0186223280616105 3305995.3178067337721586, 4066499.5436093052849174 3307244.3738269302994013, 4064267.5918585541658103 3308612.1156188808381557, 4062110.2809946816414595 3310094.7942974045872688, 4060033.5240568136796355 3311688.3459433037787676, 4058043.0132913626730442 3313388.4027422824874520, 4056144.2045499780215323 3315190.3049568249844015, 4054342.3023354355245829 3317089.1136982091702521, 4052642.2455364568158984 3319079.6244636601768434, 4051048.6938905580900609 3321156.3814015286043286, 4049566.0152120343409479 3323313.6922654006630182, 4048198.2734200833365321 3325545.6440161522477865, 4046949.2173998868092895 3327846.1190291745588183, 4045822.2707271752879024 3330208.8118623620830476, 4044820.5222844453528523 3332627.2465388872660697, 4043946.7177945477887988 3335094.7942974045872688, 4043559.9463409590534866 3336400.5128421927802265, 4044044.7568011740222573 3335031.4503365186974406, 4045046.5052439039573073 3332613.0156599935144186, 4046173.4519166154786944 3330250.3228268059901893, 4047422.5079368120059371 3327949.8478137836791575, 4048790.2497287630103528 3325717.8960630320943892, 4049331.1645958102308214 3324930.8600477352738380)), POLYGON ((4043559.9463409590534866 3336400.5128421927802265, 4043203.2522948519326746 3337604.6917610261589289, 4043559.9463409585878253 3336400.5128421946428716, 4043559.9463409590534866 3336400.5128421927802265)), POLYGON ((4048198.2734200833365321 3375545.6440161522477865, 4049566.0152120338752866 3377777.5957669033668935, 4048532.4282185747288167 3376090.9350878326222301, 4048532.4282185737974942 3376090.9350878312252462, 4048198.2734200833365321 3375545.6440161522477865)), POLYGON ((4048677.9301270744763315 3375996.8418176323175430, 4047428.8741068779490888 3373696.3668046095408499, 4046301.9274341664277017 3371333.6739714220166206, 4047428.8741068779490888 3373696.3668046100065112, 4048677.9301270744763315 3375996.8418176323175430)), POLYGON ((4045300.1789914364926517 3368915.2392948972992599, 4044426.3745015384629369 3366447.6915363795123994, 4043682.9090018430724740 3363937.7940727584064007, 4044426.3745015379972756 3366447.6915363795123994, 4045300.1789914364926517 3368915.2392948972992599)), POLYGON ((4042052.5128941014409065 3353674.4936534985899925, 4042061.0430040201172233 3354000.2454089340753853, 4042266.4249733351171017 3356609.8707601693458855, 4042061.0430040205828846 3354000.2454089340753853, 4042052.5128941014409065 3353674.4936534985899925)), POLYGON ((4046679.1870039231143892 3372483.2961305696517229, 4047442.1935323304496706 3374082.9725837642326951, 4048532.4282185714691877 3376090.9350878265686333, 4048532.4282185742631555 3376090.9350878312252462, 4047442.1935323309153318 3374082.9725837642326951, 4046679.1870039231143892 3372483.2961305696517229)))
# </A>
# Error in CPL_nary_intersection(x) : GEOS exception
# Calls: %>% ... st_intersection -> st_intersection.sfc -> CPL_nary_intersection
# Execution halted

With setting precision to 1, I see

library(sf)
# Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(dplyr)

# Attaching package: ‘dplyr’

# The following objects are masked from ‘package:stats’:

#     filter, lag

# The following objects are masked from ‘package:base’:

#     intersect, setdiff, setequal, union

df <- structure(list(lat = c(53.218461547511, 53.2354543814792, 53.2219368310548,
53.2142084420887, 53.2051969857393), lon = c(6.57022923569802,
6.55005158685787, 6.57015231400986, 6.5633788874466, 6.5692156138519
)), row.names = c(NA, -5L), class = c("tbl_df", "tbl", "data.frame"
))

#Then I create a circle around each point:

circle_df <- df %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
st_transform(3035) %>%
st_buffer(dist = units::set_units(50, "kilometers"))

#Applying st_intersection() gives an error:

y <- circle_df %>% st_set_precision(1) %>% st_intersection()
# geometry errors: 1
# Warning message:
# In CPL_nary_intersection(x) : one or more notices ignored

@MHaringa
Copy link
Author

MHaringa commented May 3, 2021

Because GEOS runs into an exception. You created a very difficult problem for it, with data like this:
x

Thanks for your quick response. Can you elaborate a bit more why my problem is very difficult for GEOS? Setting precision to 1 gives no errors indeed. Can I also get around the error by removing points that are very close to other points?

@dbaston
Copy link
Contributor

dbaston commented May 3, 2021

It should work without setting precision, but to find out why it's not we'd need to pull out the inputs to GEOS that are failing.

@edzer
Copy link
Member

edzer commented May 3, 2021

As WKT? :

options(digits = 20)
circle_df %>% st_geometry() %>% st_as_text()

@dbaston
Copy link
Contributor

dbaston commented May 3, 2021

Well, I mean the geometries to the GEOS function, not the sf function. So you (or me, later on) probably needs to patch the code to dump those to Rcout ?

@dbaston
Copy link
Contributor

dbaston commented May 4, 2021

It looks like what's happening here is that one of the processing steps produces a mixed dimension output (a GeometryCollection with LineString and Polygon components.) That causes a subsequent call to GEOSDifference to fail ("Overlay input is mixed-dimension.") I guess GEOS needs some way to only retain the highest-dimension outputs in a case like this. (cc @dr-jts)

@dr-jts
Copy link

dr-jts commented May 4, 2021

Yes, when overlay operations return mixed results it causes problems for chained operations as in this situation.

This was a long-standing problem with the original overlay algorithm (from my perspective). I was hoping to eliminate it in OverlayNG by returning only homogeneous resuts where applicable. But this was not popular, due to the semantic change. So the old behaviour was maintained as the default. But there is an OverlayNG strict mode, which does have this behaviour. However, it isn't surfaced up to OverlayNGRobust, which it really needs to be to be useful.

I see that GEOS did port the strict mode flag, so it's at least partway there too.

In the meantime, the best thing to do is to write a utility function which keeps only the highest dimension elements in the overlay result, for certain operations (intersection and difference).

@dbaston
Copy link
Contributor

dbaston commented May 4, 2021

Ticketed at https://trac.osgeo.org/geos/ticket/1117. Unfortunately I don't think there's anything that can be done for the 3.9 series. And adding a function for 3.10 will make for some awkward #ifdef-ing.

@edzer
Copy link
Member

edzer commented May 4, 2021

Thanks, @dbaston ! If this is not going to happen any time soon in GEOS, we could also follow the suggestion of Martin @dr-jts and get rid of lower dimensional geometries in a collection in the loop(s) where they're formed, in this package.

@vincentmaarek
Copy link

Hi everyone,
Does anyone have an update about this issue ?
I face it quite often and I have still not been able to find a proper workaround.

@edzer, the idea of getting rid of lower dimensional geometries in the loop would be ideal in my opinion. Not sure if this is planned in the near future, but if not, is there any possibility to share the underlying algorithm of the st_intersection() function, in order for me to write some R code that would reproduce it and that would remove lower dimensional geometries ?

Not sure to know where I can find out the source of the algorithm.
Thanks for your help.

Kind regards,
Vincent

@michal-kinel
Copy link

Using branch geos_precision, without setting precision I see errors from GEOS 3.9.0:

With setting precision to 1, I see

library(sf)
# Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(dplyr)

# Attaching package: ‘dplyr’

# The following objects are masked from ‘package:stats’:

#     filter, lag

# The following objects are masked from ‘package:base’:

#     intersect, setdiff, setequal, union

df <- structure(list(lat = c(53.218461547511, 53.2354543814792, 53.2219368310548,
53.2142084420887, 53.2051969857393), lon = c(6.57022923569802,
6.55005158685787, 6.57015231400986, 6.5633788874466, 6.5692156138519
)), row.names = c(NA, -5L), class = c("tbl_df", "tbl", "data.frame"
))

#Then I create a circle around each point:

circle_df <- df %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
st_transform(3035) %>%
st_buffer(dist = units::set_units(50, "kilometers"))

#Applying st_intersection() gives an error:

y <- circle_df %>% st_set_precision(1) %>% st_intersection()
# geometry errors: 1
# Warning message:
# In CPL_nary_intersection(x) : one or more notices ignored

in this case this worked
y <- circle_df %>% lwgeom::st_snap_to_grid(50) %>% st_set_precision(1) %>% st_intersection()

@edzer
Copy link
Member

edzer commented Nov 11, 2021

Thanks; I don't think that branch exists anymore.

@edzer
Copy link
Member

edzer commented Apr 7, 2023

ad34f46 uses the GEOS overlayNG routines, ending in ...Prec_r() for (GEOS) Intersection, Difference, Union and SymDifference; see #2143

@jepol77
Copy link

jepol77 commented Oct 12, 2023

This should work:

st_stack <- function(x) {
require(dplyr)
require(sf)
sf_use_s2(F)

#Function does not work when there are duplicate geometries
if(length(unique(x$geometry)) != nrow(x))
stop("Error: Duplicate geometries in input data")

#Faster st_erase function
st_erase_big = function(x, y){
x$ref <- 1:nrow(x)
x_inter = x[apply(st_intersects(x,y), 1, any),] #x in y
y_inter = y[apply(st_intersects(y,x), 1, any),] #y in x
difference <- st_difference(x_inter, st_union(y_inter)) #x_inter not in y_inter
nointer_x = x[!x$ref %in% x_inter$ref,] # x not intersecting
dplyr::bind_rows(difference,nointer_x)
}

#Progress bar
progress <- function (x, max = 100) {
percent <- x / max * 100
cat(sprintf('\r[%-50s] %d%%',
paste(rep('=', percent / 2), collapse = ''),
floor(percent)))
if (x == max)
cat('\n')
}

a <- st_as_sf(x[1,"geometry"])
a$n_overlaps <- 1

d <- x[,"geometry"]
e <- x[,"geometry"]
e$n <- 1
i <- 2

for(i in 2:nrow(x)){
progress(match(i, 2:nrow(x)),length(2:nrow(x)))
b <- st_intersection(a, e[i,])

#If no overlaps - just bind them together
if(nrow(b) == 0){
a <- plyr::rbind.fill(list(a, d[i,])) %>%
mutate(n_overlaps= ifelse(is.na(n_overlaps), 1, n_overlaps)) %>%
st_as_sf() %>%
st_zm()

next
}

b$n_overlaps <- b$n_overlaps+b$n
b$n <- NULL

#Get geometry type
ty <- st_geometry_type(a, by_geometry = TRUE)

#if some are a collection, only take polygons from those.
if(!all(unique(ty) %in% c("POLYGON", "MULTIPOLYGON"))){
zz <- a[ty == "GEOMETRYCOLLECTION",] %>%
st_collection_extract("POLYGON")
a <- rbind(a[ty != "GEOMETRYCOLLECTION",], zz)
}

#Remove areas from new polygon
zz <- st_erase_big(d[i,], a[,"geometry"]) %>%
mutate(ref=NULL)

#Remove areas from old polygons
a <- st_erase_big(a, d[i,])%>%
mutate(ref=NULL)

#Bind old polygons, new polygon and intersected area
a <- plyr::rbind.fill(list(a, b, zz)) %>%
mutate(n_overlaps= ifelse(is.na(n_overlaps), 1, n_overlaps)) %>%
st_as_sf() %>%
st_zm()
}
return(a)
}

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

7 participants