Почему sf :: st_transform возвращает пустую геометрию при применении к df со столбцом sf (POINT)? - PullRequest
Новая
       103

Почему sf :: st_transform возвращает пустую геометрию при применении к df со столбцом sf (POINT)?

0 голосов
/ 26 апреля 2020

Я пытаюсь st_transform sf-точек (CRS: 4326, в градусах) в кадре данных в локальную метри c проекцию (EPSG: 23867, DGN95 / UTM зона 47N).

Итак, я есть датафреймы "waypoints" со всеми моими атрибутами, а также lon / lat как столбцы. Я создаю точечную геометрию sf в кадре данных с "lat" и "lon" в качестве входных данных (waypoints_sf). (дополнительно я сохраняю координаты как Север и Восток в двух новых столбцах) 

waypoints_sf <-  st_as_sf(waypoints,
                      coords = c("lat", "lon"),
                      crs = 4326)
coordinates_tmp <- st_coordinates(waypoints_sf)
colnames(coordinates_tmp) <- c("E","N") 
waypoints_sf <- cbind(waypoints_sf,coordinates_tmp) # save coordinates as column E and N

class(waypoints_sf)
[1] "sf"         "grouped_df" "tbl_df"     "tbl"        "data.frame"

class(waypoints_sf$geometry)
[1] "sfc_POINT" "sfc"

st_crs(waypoints_sf$geometry)
[1] Coordinate Reference System:
    EPSG: 4326 
    proj4string: "+proj=longlat +datum=WGS84 +no_defs"

Теперь, когда я пытаюсь st_transform этот фрейм данных, который не выдает ошибку или предупреждение, но точки sf являются пустыми геометриями c (NaN, NaN) что бы я ни пытался. 

 waypoints_sf <-  st_transform(waypoints_sf, 23867)
 waypoints_sf$geometry[1]
 [1] Geometry set for 1 feature  (with 1 geometry empty)
     geometry type:  POINT
     dimension:      XY
     bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
     epsg (SRID):    23867
     proj4string:    +proj=utm +zone=47 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
     POINT EMPTY

Итак, что я делаю не так, может быть, я действительно скучаю по простому, я действительно не знаю, почему это не работает? Какой будет простой обходной путь, если он не работает с st_transform, возможно, с использованием N и E raw?

Обходной путь, который на данный момент работает для меня ("waypoints" - это простой df): 

waypoints_sf <- waypoints
waypoints_sf <- waypoints_sf %>% mutate(longitude=lon,latitude=lat) # save if later needed
coordinates(waypoints_sf) <- c("lon","lat")
proj4string(waypoints_sf) <- CRS("+proj=longlat +datum=WGS84") # set crs
waypoints_sf <- spTransform(test, CRS("+init=epsg:23867")) # transform
waypoints_sf <- st_as_sf(waypoints_sf,crs = 23867)

1 Ответ

2 голосов
/ 26 апреля 2020

Я подозреваю, что у вас поменялись широта и долгота, что приводит к недопустимой комбинации в контексте EPSG: 23867

Рассмотрите этот код (обратите внимание, что я поменял ваши координаты широты!) 

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

waypoints_sf <-  readr::read_csv("sample.csv") %>% 
  st_as_sf(coords = c("lon", "lat"), # note the order!
           crs = 4326)

# this works...
waypoints_trans <-  st_transform(waypoints_sf, crs = 23867)

# sanity check - a leaflet plot of the original sf data frame
leaflet(data = waypoints_sf) %>% 
  addTiles() %>% 
  addCircleMarkers(fillColor = "red",
                   stroke = FALSE)

А в качестве проверки работоспособности и окончательного подтверждения карта в листовке; это выглядит хорошо? Да, тогда ваши координаты были перевернуты.

enter image description here

...