英文:
sf::st_intersection not implement all points
问题
我有两个sf文件：x 和 Riv
我想要获取 'x' 线段与 'Riv' 线段的交集。为此，我使用了函数 `sf::st_intersection(x, Riv)`。
然而，输出结果并不包括所有 'x' 线段与 'Riv' 的交集。
(只有 27 个交集中的 44 个)
sf文件 'x' 是通过在Riv和一个早期的sf文件（在下面的图中以绿色显示）之间应用 `st_nearest_points()` 函数得到的。然后，所有 'x' 线段应该与 'Riv' 有交集。
我想找到一种方法，将所有缺失的 'x' 线段包括在输出结果中（在图上以红色显示），而不仅仅是几个被选中的（在图上以蓝色显示）。
x <- structure(list(structure(c(163340.962516503, 163339.007569445,
86055.2969082572, 86057.5225402927), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), ...), n_empty = 0L, precision = 0, crs = structure(list(
input = "BD72 / Belgian Lambert 72", wkt = "PROJCRS["BD72 / Belgian Lambert 72",\n BASEGEOGCRS["BD72",\n DATUM["Reseau National Belge 1972",\n ELLIPSOID["International 1924",6378388,297,\n LENGTHUNIT["metre",1]]],\n PRIMEM["Greenwich",0,\n ANGLEUNIT["degree",0.0174532925199433]],\n ID["EPSG",4313]],\n CONVERSION["Belgian Lambert 72",\n METHOD["Lambert Conic Conformal (2SP)",\n ID["EPSG",9802]],\n PARAMETER["Latitude of false origin",90,\n ANGLEUNIT["degree",0.0174532925199433],\n ID["EPSG",8821]],\n PARAMETER["Longitude of false origin",4.36748666666667,\n ANGLEUNIT["degree",0.0174532925199433],\n ID["EPSG",8822]],\n PARAMETER["Latitude of 1st standard parallel",51.1666672333333,\n ANGLEUNIT["degree",0.0174532925199433],\n ID["EPSG",8823]],\n PARAMETER["Latitude of 2nd standard parallel",49.8333339,\n ANGLEUNIT["degree",0.0174532925199433],\n ID["EPSG",8824]],\n PARAMETER["Easting at false origin",150000.013,\n LENGTHUNIT["metre",1],\n ID["EPSG",8826]],\n PARAMETER["Northing at false origin",5400088.438,\n LENGTHUNIT["metre",1],\n ID["EPSG",8827]]],\n CS[Cartesian,2],\n AXIS["easting (X)",east,\n ORDER[1],\n LENGTHUNIT["metre",1]],\n AXIS["northing (Y)",north,\n ORDER[2],\n LENGTHUNIT["metre",1]],\n USAGE[\n SCOPE["Engineering survey, topographic mapping."],\n AREA["Belgium - onshore."],\n BBOX[49.5,2.5,51.51,6.4]],\n ID["EPSG",31370]]"), class = "crs"), class = c("sfc_LINESTRING",
"sfc"), bbox = structure(c(xmin = 163319.91809282, ymin = 86046.2351102559,
xmax = 163632.39098375, ymax = 86190.762703538), class = "bbox"))
Riv <- structure(list(id = 1, Type = NA_character_, geometry = structure(list(
structure(c(163293.837499969, 163298.070833307, 163303.18611109,
...), dim = c(54L, 2L), class = c("
英文:
I have two sf files: x & Riv
I would like to obtain the intersection of 'x' polylines with 'Riv' Linestring. For that, I used the function sf::st_intersection(x,Riv).
Nevertheless, the output does not include all the intersection with all 'x' linestrings.
(Only 27 on 44)
The sf file 'x' come from the application of the st_nearest_points()between Riv and a former sf file (in green in the figure on below). Then all 'x' lines should have an intersection with 'Riv'.
I would like to find a way to include all the missing x lines in the outputs (in red in the on figure) and not only the few selected (in blue on the figure).
x <- structure(list(structure(c(163340.962516503, 163339.007569445,
86055.2969082572, 86057.5225402927), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163344.400724905, 163343.429912257,
86058.3169561776, 86059.4221890382), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163343.975166494, 163343.159782333,
86057.9431548709, 86058.8714383775), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163340.371658435, 163338.188246187,
86054.7779113055, 86057.2636421723), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163358.985477678, 163358.405840557,
86057.1043793742, 86056.243204223), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163337.034088493, 163336.815430932,
86052.5278108928, 86052.9286830882), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163343.251330057, 163341.479437086,
86057.3073525952, 86059.3245845931), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163340.38060113, 163339.576337948,
86054.785766375, 86055.7013890743), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163342.6099935, 163341.283701168,
86056.7440164302, 86058.2539492389), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163336.398516082, 163335.875658191,
86052.1811350321, 86053.1397078326), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163359.765722677, 163357.848655319,
86056.5823423227, 86053.7067412864), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163341.887132772, 163341.429950661,
86056.109071196, 86056.6295554452), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163341.612846091, 163341.473399649,
86055.8681437061, 86056.0268981168), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163337.333727026, 163336.984452273,
86052.6912500924, 86053.3315871395), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163347.378671305, 163347.203945911,
86060.0021536016, 86060.4972088849), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163330.103100902, 163330.035532865,
86049.5971912096, 86050.0279374449), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163351.681033356, 163351.828097384,
86060.6313092021, 86061.3265209692), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163358.810584269, 163359.400106365,
86057.2220960919, 86058.0979574919), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163320.051020384, 163319.91809282,
86049.0132598105, 86048.122645136), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163322.057491585, 163321.958762763,
86048.7137864969, 86048.0523033906), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163321.663260317, 163321.531477757,
86048.7726269846, 86047.8896838333), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163321.720228212, 163321.388596349,
86048.7641243137, 86046.5421908274), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163321.179157688, 163321.200699036,
86048.8448811083, 86048.9892081348), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163355.863493616, 163354.256530111,
86059.2057148, 86056.8182261633), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163321.776604765, 163321.528371423,
86048.7557099029, 86047.0925465152), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163320.763208813, 163320.771370347,
86048.90696303, 86048.9616453042), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163358.156758331, 163357.230257917,
86057.662171242, 86056.2856563414), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163329.134219573, 163328.981016192,
86049.4452098247, 86050.4218813749), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163321.814856331, 163321.834515089,
86048.7500007138, 86048.8817143925), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163320.297973082, 163320.356287921,
86048.9764011989, 86049.3671106212), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163321.747157576, 163321.715059441,
86048.7601050058, 86048.545047503), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163356.166413244, 163354.726073394,
86059.001826589, 86056.8618930979), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163363.182741817, 163363.228177877,
86054.4241853348, 86055.0754355295), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163514.131668895, 163514.116761467,
86092.5347591268, 86092.984963432), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163536.571089144, 163536.637359745,
86093.7360065682, 86092.5002547856), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163423.973839794, 163424.105831473,
86059.3596477238, 86058.9129066588), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163632.368401207, 163632.39098375,
86190.762703538, 86190.4820347959), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163444.094941805, 163443.50918856,
86060.4176706545, 86062.4491127571), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163360.049154527, 163359.108284222,
86056.3933877558, 86054.982082298), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163356.474576667, 163355.752668853,
86058.7944089005, 86057.721860148), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163362.583680115, 163361.324256662,
86054.7037040303, 86052.8145688502), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163413.172675529, 163413.167783795,
86057.2378796361, 86058.0401240094), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163320.98844555, 163320.710672094,
86048.8733456066, 86047.0122634526), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(163321.10674043, 163320.715609176,
86048.8556896544, 86046.2351102559), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg"))), n_empty = 0L, precision = 0, crs = structure(list(
input = "BD72 / Belgian Lambert 72", wkt = "PROJCRS[\"BD72 / Belgian Lambert 72\",\n BASEGEOGCRS[\"BD72\",\n DATUM[\"Reseau National Belge 1972\",\n ELLIPSOID[\"International 1924\",6378388,297,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4313]],\n CONVERSION[\"Belgian Lambert 72\",\n METHOD[\"Lambert Conic Conformal (2SP)\",\n ID[\"EPSG\",9802]],\n PARAMETER[\"Latitude of false origin\",90,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8821]],\n PARAMETER[\"Longitude of false origin\",4.36748666666667,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8822]],\n PARAMETER[\"Latitude of 1st standard parallel\",51.1666672333333,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8823]],\n PARAMETER[\"Latitude of 2nd standard parallel\",49.8333339,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8824]],\n PARAMETER[\"Easting at false origin\",150000.013,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8826]],\n PARAMETER[\"Northing at false origin\",5400088.438,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8827]]],\n CS[Cartesian,2],\n AXIS[\"easting (X)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"northing (Y)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"Belgium - onshore.\"],\n BBOX[49.5,2.5,51.51,6.4]],\n ID[\"EPSG\",31370]]"), class = "crs"), class = c("sfc_LINESTRING",
"sfc"), bbox = structure(c(xmin = 163319.91809282, ymin = 86046.2351102559,
xmax = 163632.39098375, ymax = 86190.762703538), class = "bbox"))
Riv <- Riv <- structure(list(id = 1, Type = NA_character_, geometry = structure(list(
structure(c(163293.837499969, 163298.070833307, 163303.18611109,
163310.77083332, 163317.38541666, 163323.294444444, 163332.290277786,
163339.081250015, 163345.607638911, 163350.105555582, 163354.691666698,
163359.277777814, 163362.981944484, 163366.774305599, 163371.007638937,
163373.477083384, 163416.868750095, 163424.629861214, 163427.716666772,
163443.238889011, 163457.61458347, 163468.021527925, 163473.136805708,
163479.839583492, 163489.71736128, 163500.918055736, 163514.235416861,
163542.193055778, 163550.659722453, 163557.450694683, 163561.772222465,
163570.503472474, 163574.031250255, 163573.590278032, 163570.679861363,
163566.181944692, 163560.802083575, 163561.419444687, 163563.359722466,
163568.29861136, 163569.533333584, 163573.325694699, 163575.706944701,
163583.820833598, 163593.610416942, 163595.9034725, 163599.960416948,
163605.957639177, 163622.802778083, 163638.148611432, 163646.174305884,
163659.668055898, 163670.251389242, 163690.536111485, 86050.9986111142,
86054.438194451, 86056.4666666753, 86053.4680555611, 86049.4111111126,
86048.5291666672, 86049.9402777798, 86053.6444444502, 86059.377083345,
86060.9645833466, 86059.9944444567, 86056.907638898, 86054.438194451,
86054.1736111174, 86055.761111119, 86056.9958333425, 86057.2604166761,
86059.553472234, 86060.6118055684, 86060.1708333457, 86064.3159722389,
86071.4597222462, 86078.1625000308, 86087.2465278178, 86090.6861111547,
86092.0972222672, 86092.5381944899, 86094.037500047, 86097.2125000502,
86099.5937500526, 86104.444444502, 86116.8798611814, 86127.375000081,
86133.9013889765, 86141.3979167619, 86146.7777778785, 86151.7166667725,
86155.9500001101, 86165.3868056753, 86181.7909723587, 86187.7000001425,
86196.8722223741, 86201.2819446008, 86203.8395834923, 86203.8395834923,
86200.4000001555, 86195.372916817, 86192.5506945919, 86189.9930557004,
86191.2277779239, 86193.4326390372, 86194.1381945935, 86193.7854168154,
86186.5534723636), dim = c(54L, 2L), class = c("XY", "LINESTRING",
"sfg"))), n_empty = 0L, crs = structure(list(input = "BD72 / Belgian Lambert 72",
wkt = "PROJCRS[\"BD72 / Belgian Lambert 72\",\n BASEGEOGCRS[\"BD72\",\n DATUM[\"Reseau National Belge 1972\",\n ELLIPSOID[\"International 1924\",6378388,297,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4313]],\n CONVERSION[\"Belgian Lambert 72\",\n METHOD[\"Lambert Conic Conformal (2SP)\",\n ID[\"EPSG\",9802]],\n PARAMETER[\"Latitude of false origin\",90,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8821]],\n PARAMETER[\"Longitude of false origin\",4.36748666666667,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8822]],\n PARAMETER[\"Latitude of 1st standard parallel\",51.1666672333333,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8823]],\n PARAMETER[\"Latitude of 2nd standard parallel\",49.8333339,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8824]],\n PARAMETER[\"Easting at false origin\",150000.013,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8826]],\n PARAMETER[\"Northing at false origin\",5400088.438,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8827]]],\n CS[Cartesian,2],\n AXIS[\"easting (X)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"northing (Y)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"Belgium - onshore.\"],\n BBOX[49.5,2.5,51.51,6.4]],\n ID[\"EPSG\",31370]]"), class = "crs"), class = c("sfc_LINESTRING",
"sfc"), precision = 0, bbox = structure(c(xmin = 163293.837499969,
ymin = 86048.5291666672, xmax = 163690.536111485, ymax = 86203.8395834923
), class = "bbox"))), row.names = c(NA, -1L), class = c("sf",
"data.frame"), sf_column = "geometry", agr = structure(c(id = NA_integer_,
Type = NA_integer_), class = "factor", levels = c("constant",
"aggregate", "identity")))
EDIT: The original sf object, used for st_nearest_point()
ori <- structure(list(N = c(1, 2, 3, 10, 12, 20, 23, 25, 26, 30, 34,
46, 51, 61, 62, 68, 69, 83, 2, 4, 6, 7, 8, 14, 17, 19, 21, 24,
25, 26, 30, 32, 36, 37, 40, 44, 55, 57, 58, 59, 74, 3, 31, 32
), geometry = structure(list(structure(c(163339.007569445, 86057.5225402927
), class = c("XY", "POINT", "sfg")), structure(c(163343.429912257,
86059.4221890382), class = c("XY", "POINT", "sfg")), structure(c(163343.159782333,
86058.8714383775), class = c("XY", "POINT", "sfg")), structure(c(163338.188246187,
86057.2636421723), class = c("XY", "POINT", "sfg")), structure(c(163358.405840557,
86056.243204223), class = c("XY", "POINT", "sfg")), structure(c(163336.815430932,
86052.9286830882), class = c("XY", "POINT", "sfg")), structure(c(163341.479437086,
86059.3245845931), class = c("XY", "POINT", "sfg")), structure(c(163339.576337948,
86055.7013890743), class = c("XY", "POINT", "sfg")), structure(c(163341.283701168,
86058.2539492389), class = c("XY", "POINT", "sfg")), structure(c(163335.875658191,
86053.1397078326), class = c("XY", "POINT", "sfg")), structure(c(163357.848655319,
86053.7067412864), class = c("XY", "POINT", "sfg")), structure(c(163341.429950661,
86056.6295554452), class = c("XY", "POINT", "sfg")), structure(c(163341.473399649,
86056.0268981168), class = c("XY", "POINT", "sfg")), structure(c(163336.984452273,
86053.3315871395), class = c("XY", "POINT", "sfg")), structure(c(163347.203945911,
86060.4972088849), class = c("XY", "POINT", "sfg")), structure(c(163330.035532865,
86050.0279374449), class = c("XY", "POINT", "sfg")), structure(c(163351.828097384,
86061.3265209692), class = c("XY", "POINT", "sfg")), structure(c(163359.400106365,
86058.0979574919), class = c("XY", "POINT", "sfg")), structure(c(163319.91809282,
86048.122645136), class = c("XY", "POINT", "sfg")), structure(c(163321.958762763,
86048.0523033906), class = c("XY", "POINT", "sfg")), structure(c(163321.531477757,
86047.8896838333), class = c("XY", "POINT", "sfg")), structure(c(163321.388596349,
86046.5421908274), class = c("XY", "POINT", "sfg")), structure(c(163321.200699036,
86048.9892081348), class = c("XY", "POINT", "sfg")), structure(c(163354.256530111,
86056.8182261633), class = c("XY", "POINT", "sfg")), structure(c(163321.528371423,
86047.0925465152), class = c("XY", "POINT", "sfg")), structure(c(163320.771370347,
86048.9616453042), class = c("XY", "POINT", "sfg")), structure(c(163357.230257917,
86056.2856563414), class = c("XY", "POINT", "sfg")), structure(c(163328.981016192,
86050.4218813749), class = c("XY", "POINT", "sfg")), structure(c(163321.834515089,
86048.8817143925), class = c("XY", "POINT", "sfg")), structure(c(163320.356287921,
86049.3671106212), class = c("XY", "POINT", "sfg")), structure(c(163321.715059441,
86048.545047503), class = c("XY", "POINT", "sfg")), structure(c(163354.726073394,
86056.8618930979), class = c("XY", "POINT", "sfg")), structure(c(163363.228177877,
86055.0754355295), class = c("XY", "POINT", "sfg")), structure(c(163514.116761467,
86092.984963432), class = c("XY", "POINT", "sfg")), structure(c(163536.637359745,
86092.5002547856), class = c("XY", "POINT", "sfg")), structure(c(163424.105831473,
86058.9129066588), class = c("XY", "POINT", "sfg")), structure(c(163632.39098375,
86190.4820347959), class = c("XY", "POINT", "sfg")), structure(c(163443.50918856,
86062.4491127571), class = c("XY", "POINT", "sfg")), structure(c(163359.108284222,
86054.982082298), class = c("XY", "POINT", "sfg")), structure(c(163355.752668853,
86057.721860148), class = c("XY", "POINT", "sfg")), structure(c(163361.324256662,
86052.8145688502), class = c("XY", "POINT", "sfg")), structure(c(163413.167783795,
86058.0401240094), class = c("XY", "POINT", "sfg")), structure(c(163320.710672094,
86047.0122634526), class = c("XY", "POINT", "sfg")), structure(c(163320.715609176,
86046.2351102559), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",
"sfc"), precision = 0, bbox = structure(c(xmin = 163319.91809282,
ymin = 86046.2351102559, xmax = 163632.39098375, ymax = 86190.4820347959
), class = "bbox"), crs = structure(list(input = "EPSG:31370",
wkt = "PROJCRS[\"BD72 / Belgian Lambert 72\",\n BASEGEOGCRS[\"BD72\",\n DATUM[\"Reseau National Belge 1972\",\n ELLIPSOID[\"International 1924\",6378388,297,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4313]],\n CONVERSION[\"Belgian Lambert 72\",\n METHOD[\"Lambert Conic Conformal (2SP)\",\n ID[\"EPSG\",9802]],\n PARAMETER[\"Latitude of false origin\",90,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8821]],\n PARAMETER[\"Longitude of false origin\",4.36748666666667,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8822]],\n PARAMETER[\"Latitude of 1st standard parallel\",51.1666672333333,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8823]],\n PARAMETER[\"Latitude of 2nd standard parallel\",49.8333339,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8824]],\n PARAMETER[\"Easting at false origin\",150000.013,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8826]],\n PARAMETER[\"Northing at false origin\",5400088.438,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8827]]],\n CS[Cartesian,2],\n AXIS[\"easting (X)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"northing (Y)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"Belgium - onshore.\"],\n BBOX[49.5,2.5,51.51,6.4]],\n ID[\"EPSG\",31370]]"), class = "crs"), n_empty = 0L)), row.names = c(NA,
44L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(N = NA_integer_), levels = c("constant",
"aggregate", "identity"), class = "factor"))
答案1
得分: 2
以下是翻译好的内容：
There are missing intersections, so those failing x lines are still slightly off, with gap distances in ~ 1e-10 range. Precision and floating point representation can hit hard ... experimenting with precision and rounding could work, but let's test sfnetworks. And let's use points from ori, the original sf object.
注意，这将改变几何形状 - Riv 在每个混合/交叉点处被分割，然后联合在一起，结果形状将不完全相同于原始的。而且我无法猜测它在性能方面是否能扩展。
library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(sfnetworks)
library(ggplot2)
library(patchwork)
library(dplyr)
# 检查 `x` 和 Riv 之间的距离，单位为米
# 一切 > 0 将不相交
st_distance(x, Riv) |> c()
#> Units: [m]
#> [1] 5.127934e-11 0.000000e+00 1.560592e-10 0.000000e+00 0.000000e+00
#> ...
#> [41] 5.415751e-11 3.582375e-11 1.468602e-11 0.000000e+00
# sfnetworks 提供方便的 st_network_blend()
net <- as_sfnetwork(Riv) |>
st_network_blend(ori)
autoplot(net) + theme_bw()
# 从网络中提取节点（混合点），过滤掉 Riv 的端点
n_points <- net |>
activate(nodes) |>
filter(!is.na(N)) |>
st_as_sf()
# 从网络中提取边（分割的 Riv，线段），联合为多线
n_riv <- net |>
activate(edges) |>
st_as_sf() |>
st_union()
# 检查距离
st_distance(n_points, n_riv) |> c()
#> Units: [m]
#> [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [39] 0 0 0 0 0 0
# 检查相交点
st_intersects(n_points, n_riv) |> t() |> unlist()
#> [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
#> [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
ggplot() +
geom_sf(data = n_riv, alpha = .5) +
geom_sf(data = n_points, aes(color = "n_points"), shape = 4) +
geom_sf(data = ori, aes(color = "ori"), shape = 4) +
geom_sf(data = x, aes(color = "x")) +
theme(legend.position = "none") +
theme_bw()
autoplot(net):
<!-- -->
结果几何形状以及参考的 x 线：
<!-- -->
<!-- -->
<!-- -->
英文:
There are missing intersections, so those failing x lines are still slightly off, with gap distances in ~ 1e-10 range. Precision and floating point representation can hit hard ... experimenting with precision and rounding could work, but let's test sfnetworks. And let's use points from ori, the original sf object.
Note that this will alter geometries -- Riv is split at each blend / intersection point and later unioned back together, resulting shape will not be identical to the original. And I can't even guess if/how it scales performance-wise.
library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(sfnetworks)
library(ggplot2)
library(patchwork)
library(dplyr)
# check distances between x and Riv, meters
# everything > 0 will not intersect
st_distance(x,Riv) |> c()
#> Units: [m]
#> [1] 5.127934e-11 0.000000e+00 1.560592e-10 0.000000e+00 0.000000e+00
#> ...
#> [41] 5.415751e-11 3.582375e-11 1.468602e-11 0.000000e+00
# sfnetworks provides convenient st_network_blend()
net <- as_sfnetwork(Riv) |>
st_network_blend(ori)
autoplot(net) + theme_bw()
# extract nodes (blend points) from network,
# filter out Riv endpoints
n_points <- net |>
activate(nodes) |>
filter(!is.na(N)) |>
st_as_sf()
# extract edges (segmented Riv, lines) from network,
# union to multiline
n_riv <- net |>
activate(edges) |>
st_as_sf() |>
st_union()
# check distance
st_distance(n_points, n_riv) |> c()
#> Units: [m]
#> [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [39] 0 0 0 0 0 0
# check intersecting points
st_intersects(n_points, n_riv) |> t() |> unlist()
#> [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
#> [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
ggplot() +
geom_sf(data = n_riv, alpha = .5) +
geom_sf(data = n_points, aes(color = "n_points"), shape = 4) +
geom_sf(data = ori, aes(color = "ori"), shape = 4) +
geom_sf(data = x, aes(color = "x")) +
theme(legend.position = "none") +
theme_bw()
autoplot(net):
<!-- -->
Resulting geometries along with x lines for reference:
<!-- -->
<!-- -->
<!-- -->
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论