geojson precision

[Robinlovelace]

Thanks for an awesome package.

I've used it for work and all went well until I looked at the object that I created.

I definitely saved a polygon representing the roof space of this site plan but plot(r_new2$finished$geometry) came out like this (it should be a nice curve following the roof which was what I drew):

Any ideas why? Have others reported this - can retry if I discover what I'm doing wrong (I did some re-editing).

Here's some partitially reproducible code, which geo-referenced the raster - that was very useful!

https://github.com/Robinlovelace/geocompr/pull/77

[timelyportfolio]

That is very unfortunate. I will look into this problem. For reference, does it seem to be a precision issue? loss of points/nodes problem?

[Robinlovelace]

Trying to reproduce and save the result in a gif. Example of not reproducing it:

[Robinlovelace]

Update: I think I can reproduce the issue and suspect it's due to zooming. .gif example to follow...

[Robinlovelace]

Maybe panning and zooming in combination:

[Robinlovelace]

[Robinlovelace]

Update: I could reproduce the issue only panning at a single zoom level so now suspect it's more about panning.

[Robinlovelace]

Update: I think this issue only occurs at high zoom levels: that's why I cannot reproduce in the first .gif example. Just tried the same thing at a high zoom level and the result looks garbled.

[tim-salabim]

Yes I can confirm. Whatever is Polygons drawn beyond a zoom level of about 15 seem garbled. Lines and points seems to be equally displaced. @timelyportfolio does this confirm your suspicion of a coordinate precision problem?

[timelyportfolio]

Strange. Here was my first test. I found a random lake in Canada and traced at zoom level 14 in Chrome. I then used editFeatures to try to trace at zoom level 16. Here is the result (note I got lazy and did not complete the trace).

Test 1 in Chrome

library(sf)
library(mapview)
library(mapedit)

# at zoom level 14 create a polygon with mapedit
#   to use for testing
p1 <- structure(list(structure(list(structure(c(-76.1935, -76.1922, 
  -76.191, -76.1904, -76.1902, -76.19, -76.1891, -76.1876, -76.1864, 
  -76.1862, -76.1862, -76.1862, -76.1857, -76.1848, -76.1834, -76.1819, 
  -76.1806, -76.1783, -76.1774, -76.1764, -76.1756, -76.1745, -76.1728, 
  -76.1722, -76.1712, -76.1699, -76.1692, -76.168, -76.1662, -76.1657, 
  -76.1647, -76.1638, -76.1627, -76.1616, -76.1608, -76.1596, -76.1595, 
  -76.1595, -76.1586, -76.1573, -76.1574, -76.1568, -76.1569, -76.1569, 
  -76.156, -76.1563, -76.1566, -76.1572, -76.1583, -76.1589, -76.16, 
  -76.161, -76.162, -76.1631, -76.1639, -76.1651, -76.167, -76.1677, 
  -76.1686, -76.1691, -76.1696, -76.1696, -76.1698, -76.1705, -76.1709, 
  -76.1712, -76.1711, -76.1715, -76.1743, -76.1777, -76.1803, -76.1813, 
  -76.182, -76.1828, -76.1844, -76.186, -76.1875, -76.1891, -76.1904, 
  -76.1918, -76.1936, -76.1949, -76.1959, -76.1968, -76.1971, -76.1976, 
  -76.1975, -76.197, -76.1962, -76.1956, -76.1951, -76.1943, -76.1935, 
  45.1188, 45.1186, 45.1187, 45.1192, 45.12, 45.1205, 45.1213, 
  45.122, 45.1224, 45.1232, 45.1239, 45.1244, 45.1249, 45.1255, 
  45.1259, 45.126, 45.1257, 45.1256, 45.126, 45.1262, 45.1262, 
  45.1259, 45.1259, 45.1257, 45.1254, 45.1249, 45.1245, 45.1248, 
  45.1244, 45.1239, 45.1236, 45.1234, 45.1238, 45.1236, 45.1232, 
  45.1233, 45.1242, 45.1248, 45.1248, 45.1248, 45.1237, 45.1229, 
  45.1224, 45.121, 45.1198, 45.1187, 45.1179, 45.1174, 45.1178, 
  45.1166, 45.1156, 45.1146, 45.1141, 45.1135, 45.1133, 45.1127, 
  45.1124, 45.1128, 45.1131, 45.1138, 45.1144, 45.1156, 45.1161, 
  45.1164, 45.1151, 45.1135, 45.1124, 45.1119, 45.1116, 45.1118, 
  45.1114, 45.1121, 45.1126, 45.1128, 45.1131, 45.1131, 45.1136, 
  45.1138, 45.1141, 45.1141, 45.1141, 45.1141, 45.1144, 45.1146, 
  45.1151, 45.1157, 45.1165, 45.1171, 45.1176, 45.1181, 45.1182, 
  45.1186, 45.1188), .Dim = c(93L, 2L))), class = c("XY", "POLYGON", 
  "sfg"))), n_empty = 0L, precision = 0, class = c("sfc_POLYGON", 
  "sfc"), bbox = structure(c(-76.1976, 45.1114, -76.156, 45.1262
  ), .Names = c("xmin", "ymin", "xmax", "ymax"), class = "bbox"), crs = structure(list(
  epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), .Names = c("epsg", 
  "proj4string"), class = "crs"))

editFeatures(st_sf(p1))

Test 2 in RStudio Viewer

I suspected the issue might be RStudio Viewer, but this test seems to yield expected results.

Test 3 Zoom from 16 to 18 during draw

Then, I thought maybe the error is caused by zooming to different levels during the draw, but everything seemed fine.

Test 4 Edit Original at zoom 18

I think I finally experienced the problem when I did editFeatures(p1) and zoomed to level 18 to edit the contour. My edit was much closer to the creek or offshoot than what appears below, so I will classify this as a replication.

[tim-salabim]

Here's an addition I drew to @timelyportfolio 's example using tst = editMap(mapview(p1)). I did the drawing in RStudio viewer at zoom level 17. I tried to follow @timelyportfolio 's outline as close as possible (i.e. replicate his vertices) for a small portion outside the lake. It turns out that I apparently did a very good job as all vertices seem to be topologically identical :-)

tst = structure(list(X_leaflet_id = 87L, feature_type = "polygon", 
    geometry = structure(list(structure(list(structure(c(-76.1705, 
    -76.1698, -76.1696, -76.1696, -76.171, -76.1709, -76.1705, 
    45.1164, 45.1161, 45.1156, 45.1144, 45.1144, 45.1151, 45.1164
    ), .Dim = c(7L, 2L))), class = c("XY", "POLYGON", "sfg"))), n_empty = 0L, precision = 0, class = c("sfc_POLYGON", 
    "sfc"), bbox = structure(c(-76.171, 45.1144, -76.1696, 45.1164
    ), .Names = c("xmin", "ymin", "xmax", "ymax"), class = "bbox"), crs = structure(list(
        epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), .Names = c("epsg", 
    "proj4string"), class = "crs"))), .Names = c("X_leaflet_id", 
"feature_type", "geometry"), row.names = 1L, sf_column = "geometry", agr = structure(c(NA_integer_, 
NA_integer_), class = "factor", .Label = c("constant", "aggregate", 
"identity"), .Names = c("X_leaflet_id", "feature_type")), class = c("sf", 
"data.frame"))

mapview(p1) + tst

[Robinlovelace]

Here are some more experiments.

Reproducible example (sorry about all the pkgs!) to geocode a building just south of Willow Terrace at Zoom Level 18:

Slightly garbled output:

[Robinlovelace]

And taking more care (zoom level: 19 in Firefox):

[Robinlovelace]

The result makes me think it's a precision/snapping issue:

[Robinlovelace]

Please try to reproduce this simple example with the code below:

library(sf)
library(mapview)
library(mapedit)
library(tmaptools)
library(stplanr)
b = tmaptools::bb("University of Leeds")
p = bb2poly(b)
sp::proj4string(p) = "+init=epsg:4326"
m = mapview(p)
p = editMap(m@map)
plot(p$finished)

[timelyportfolio]

@Robinlovelace, I used your example polygon with editFeatures(st_as_sf(p)) to try a couple of things.

Loss of Points ?

I thought maybe instead of precision we might be losing some nodes. So I edited the polygon and added nodes such that they filled the entire top edge.

But this worked as expected.

Precision ?

Next I tested precision as shown below.

This test failed with correct number of nodes but incorrect precision

[timelyportfolio]

Then I did this

and got

I used recorder = TRUE and the recorded feature matches the final sf, so the loss of precision is occurring in JavaScript. I think the newest Leaflet.draw still works in the older Leaflet. I'll try later today to use the newest Leaflet.draw to see if this might have been solved.

[Robinlovelace]

That is a great final demo of what is going on - look forward to hearing of how upstream developments can feed into the solution. Any JS people we can ask to find out why this dramatic loss of precision is happening? It sounds like you're on it in any case, many thanks for reproducing and pinning down the issue so quick.

[timelyportfolio]

I'll try the upgraded Leaflet.draw and also inspect the JavaScript GeoJSON tonight to find the source of the issue. At least, I feel like now we have an easily replicable example :)

[mdsumner]

This is such a great thread, so much useful information!

[timelyportfolio]

@Robinlovelace @tim-salabim through https://github.com/codeofsumit/leaflet.pm/pull/306#issuecomment-399096183, I just discovered https://github.com/Leaflet/Leaflet/pull/5444 which should solve or partially solve this issue. I will test.

[timelyportfolio]

@tim-salabim Leaflet > 1.0 changed geojson precision to 6 digits as described above. However, after some thorough tracing, I discovered the the st_as_sfc.geolist use of jsonlite::toJSON clipped to 4 digits which is the jsonlite default. I have changed to 6 to match Leaflet. Based on tests (shown below), this seems to be good enough to zoom levels 17, 18, and 19.

@Robinlovelace If possible, testing would be extremely helpful. Thanks!!

Example

drawing

library(sf)
library(mapview)
library(mapedit)

x <- editMap(mapview())

result

mapview(x$finished)

Questions

• Is 6 enough?

• Does the precision need to be an argument with 6 as default?

[Robinlovelace]

@timelyportfolio I'm on the case :+1:

[Robinlovelace]

@timelyportfolio what do I need to install to get the latest functionality? Just tried it and the spatial resolution issue seems to remain:

[timelyportfolio]

@Robinlovelace wow, that was quick and thanks. I think you will just need the newest (CRAN) leaflet and leaflet.extras along with devtools::install_github("timelyportfolio/mapedit"). I have not merged the pull at r-spatial yet.

[tim-salabim]

@timelyportfolio thanks so much! Awesome that it turned out to be such an easy fix. I tried various zoom levels and it seems that we loose precision at zoom 26.

IIUC, leaflet has precision 6 and that cannot be increased? Have you tried this with digits = NA in toJSON as the help states

max number of decimal digits to print for numeric values. Use I() to specify significant digits. Use NA for max precision.

which should give back the precision of what has been drawn (i.e. 6 in new leaflet), right?

Re having an argument to pass precision, I am undecided. I think most people would favor higher precision, but we cannot have more than 6, however, there may be instances where people only wanna draw at 3 or 4 digit precision, e.g. to ensure that drawn features fit into existing geometries... I guess it's easy enough to include an argument, so why not give the power to the people.

Oh, and I think we should add a link to [Wikipadia](https://en.wikipedia.org/wiki/Decimal_degrees] on this topic, as it has a nice overview on the digit -> accuracy issue (including differences at different lattitudes).

[Robinlovelace]

Confirmed: it's working! I consider this issue closed now (hence closing). I think setting the default precision to 6 is fine. Not sure how many people will want less (they could always snap after the event) but the idea of adding a precision argument sounds good. I think it's worth saying how that translates to spatial scale, e.g. objects created with mapedit have a spatial resolution of 1 millionth of a degree (0.000001), which translates to roughly 10 cm at the equator (see https://en.wikipedia.org/wiki/Decimal_degrees).

[Robinlovelace]

Gratuitous gif showing it action!

[timelyportfolio]

@tim-salabim I will explore what options we have available on the leaflet js side. Then we can decide how we would like to proceed with potential precision argument. Unfortunately, adding this will likely require a change to R leaflet or leaflet-extras, which might delay rapid incoporation. I like the idea of jsonlite NA argument which would mean we don't need to chage in future if we go with an argument.

@Robinlovelace Thanks so much for bringing this up. Sorry it took so long. I'm going to reopen for just a bit longer while we iterate through various combinations.

[Robinlovelace]

OK great, many thanks for the fix, opening the issue was the easy part! To re-iterate I think a default precision of 6 dp would be fine. Should be simpler on the coding side also. But will defer to your judgement on this - seems to have worked so far to contribute key components to an amazing package!

[timelyportfolio]

@tim-salabim your NA suggestion worked great and that way if we can get a precision argument in leaflet/leaflet.extras then we will not need to change. I worked through the code in leaflet.extras and adding a precision argument is fairly straightforward but unfortunately fairly instrusive as well. I think for now we should stick with what we have. I'll start a new issue over on leaflet.extras since that decision is more out of our hands.

Also I just realized that I didn't answer your question. leaflet > 1.1 can go to higher or lower precision, but the default is now 6.

Thanks everyone!!