Потеря точности при записи строк WKT с использованием st_as_text () в пакете R sf - PullRequest
Новая
       79

Потеря точности при записи строк WKT с использованием st_as_text () в пакете R sf

0 голосов
/ 20 марта 2020

Мне нужно повысить точность рендеринга строк wkt MULTIPOLYGON из базы данных, используя пакет sf в R.

Основой c является то, что у нас есть полигоны, нарисованные в QGIS, и WKT, вычисленные с использованием 'get wkt strings 'перед отправкой в ​​базу данных postgreSQL. Они записываются в crs 3111 (Vicgrid 94, который измеряется в метрах) и преобразуются в crs 4326 (WGS84) для отображения и отображения. Для того, чтобы получить строки в базе данных, я использую функцию st_as_text (), чтобы я мог видеть строки.

Похоже, что на этом шаге я теряю некоторую точность, поскольку некоторые координаты не достигают на 5 цифр ниже десятичного значения, как вы ожидаете для большинства картографических приложений (например, 145.101 вместо 145.10112). Я попытался использовать функцию st_precision для сброса, но она всегда на 0, что означает отсутствие округления / сжатия.

В этом текстовом состоянии wkts go в базу данных и ждет, пока я использую их.

Когда дело доходит до анализа и построения карты, эти wkts извлекаются и возвращаются к объектам sf. В приведенном ниже примере я делаю прямоугольник angular в парке, следуя по прямому пути и устанавливая несколько вершин на прямой линии на северной стороне. Этот многоугольник выглядит хорошо на протяжении всего преобразования, но выглядит ужасно (вершины отклоняются от истинного) после применения st_as_text. В некоторых случаях я теряю пару метров точности на многоугольниках.

Мне нужно иметь возможность работать с полигонами в масштабе небольших средних участков земли и, следовательно, требовать большей точности.

В приведенном ниже примере описан мой процесс, начиная от ввода строки wkt и заканчивая построением хорошей и плохой версий полигонов поверх друг друга на листовой карте.

Если бы кто-то мог помочь мне найти решение проблемы потери точности, я был бы очень признателен 

library(dplyr)
library(sf)
library(leaflet)

# WKT is imported from QGIS in the crs 3111 format
polygon <- as.data.frame( as.character( "MULTIPOLYGON (((2508897.95312214968726039 2416808.05668657133355737, 2508952.18990397034212947 2416847.28620351944118738, 2509004.84697369858622551 2416885.72586442111060023, 2509044.86634669220075011 2416914.68725277204066515, 2509082.25286619924008846 2416941.54235833324491978, 2509126.68226878298446536 2416973.53152819396927953, 2509164.39789497572928667 2417000.78156177746132016, 2509199.08573966100811958 2416930.87930171331390738, 2509229.16609074175357819 2416870.19202885124832392, 2509012.21896346053108573 2416687.47199689317494631, 2508965.35417140228673816 2416735.38993034604936838, 2508914.27681376552209258 2416790.67985356086865067, 2508897.95312214968726039 2416808.05668657133355737)))"))

# Set the name of the columnn
colnames(polygon) <- 'wkt' 

# Set the wkt column to be the sf field and tell it what the crs currently is
polygon_3111 <- 
polygon %>% 
  st_as_sf(., wkt = "wkt") %>%
  sf::`st_crs<-`(3111) 

plot(polygon_3111) # so far so good

# Change the crs to a plot-friently WGS84
polygon4326 <- 
polygon_3111 %>% 
  st_transform(4326) 

plot(polygon4326) # looks ok here too

# Now exporting that sf object as a text string. <- this is where everything goes wrong
polygon_text <-  
polygon4326$wkt %>% 
  `st_precision<-`(0) %>% # had a go at setting precision to 0
  st_as_text()

### Imagine now that the wkt field sits here in the database ready to be used later


# A few weeks later I want to start using this wkt to make maps

polygon_for_plot <- 
  polygon_text %>% 
  as.data.frame() %>% 
  `colnames<-`(.,'wkt') %>% 
  st_as_sf(., wkt = "wkt" ) %>%
  sf::`st_crs<-`(3111) 

plot(polygon_for_plot) # now it looks wonky..

# Render in a Leaflet map for a closer look

leaflet() %>% 
  addPolygons(data = polygon_for_plot, 
              color = "#444444", weight = 1, smoothFactor = 0,
              opacity = 0.0, fillOpacity = 0.3,
              fillColor = "red",
              highlightOptions = highlightOptions(color = "white", weight = 2,
                                                  bringToFront = TRUE),
              group = "polygon from wkt") %>% 
  addPolygons(data = polygon4326,
              color = "white", weight = 1, smoothFactor = 0,
              opacity = 0.0, fillOpacity = 0.3,
              fillColor = "yellow",
              highlightOptions = highlightOptions(color = "white", weight = 2,
                                                  bringToFront = TRUE),
              group = "wgs 84 polygon prior to wkt") %>% 
   addProviderTiles("Esri.WorldImagery", group = 'Esri Imagery') %>% 
  addLayersControl(
    baseGroups = c("Esri.WorldImagery"),
    overlayGroups = c("wgs 84 polygon prior to wkt",
                      "polygon from wkt"),
    options = layersControlOptions(collapsed = FALSE)
  )
...