Неизвестный CRS в QGIS при проецировании в EPSG: 25833 в R - PullRequest
Новая
       11

Неизвестный CRS в QGIS при проецировании в EPSG: 25833 в R

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

Я хочу проецировать фрейм пространственных данных в EPSG 25833 в R , но QGIS, похоже, этого не знает (для воспроизводимости я использую код jazzurro, созданный в его / ее ответе на this вопрос с небольшими изменениями) 

library(rgdal)

mydf <- structure(list(longitude = c(128.6979, 153.0046, 104.3261, 124.9019, 
                                     126.7328, 153.2439, 142.8673, 152.689), latitude = c(-7.4197, 
                                                                                          -4.7089, -6.7541, 4.7817, 2.1643, -5.65, 23.3882, -5.571)), .Names = c("longitude", 
                                                                                                                                                                 "latitude"), class = "data.frame", row.names = c(NA, -8L))


### Get long and lat from your data.frame. Make sure that the order is in lon/lat.

xy <- mydf[,c(1,2)]


# Here I use the projection EPSG:25833
spdf <- SpatialPointsDataFrame(coords = xy, data = mydf,
                               proj4string = CRS("+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))


#Export as shapefile
writeOGR(spdf, "file location", "proj_test", driver="ESRI Shapefile",overwrite_layer = T) #now I write the subsetted network as a shapefile again

Теперь, когда я загружаю шейп-файл в QGIS, он не знает проекцию.

Screenshot from QGIS

Есть идеи?

1 Ответ

2 голосов
/ 02 марта 2020

При создании SpatialPointsDataFrame: 

# Wrong!
spdf <- SpatialPointsDataFrame(coords = xy, data = mydf,
                               proj4string = CRS("+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))`

Вы указываете фрейму данных, в чем находятся ваши точки, поэтому вы должны указать 4326, так как ваши исходные данные lon / lat.

Так должно быть: 

spdf <- SpatialPointsDataFrame(coords = xy, data = mydf,
                               proj4string = CRS("+proj=longlat +datum=WGS84"))

И затем вы можете преобразовать данные в другой CRS, используя spTransform: 

spTransform(spdf, CRS('+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'))

Для этих конкретных данных мы получаем ошибку, потому что одна из точек не преобразуется в целевую CRS:

Ошибка в spTransform (xSP, CRSobj, ...): сбой в точках 3 Дополнительно: предупреждающее сообщение: в spTransform (xSP, CRSobj, ...): 1 проецируемая точка (точки) не конечна

Я предпочитаю работать в sf, поэтому мы также можем сделать: 

library(sf)
sfdf <- st_as_sf(mydf, coords = c('longitude', 'latitude'), crs=4326, remove=F)
sfdf_25833 <- sfdf %>% st_transform(25833)

sfdf_25833
#> Simple feature collection with 8 features and 2 fields (with 1 geometry empty)
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 5589731 ymin: -19294970 xmax: 11478870 ymax: 19337710
#> epsg (SRID):    25833
#> proj4string:    +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
#>   longitude latitude                   geometry
#> 1  128.6979  -7.4197 POINT (10198485 -17980649)
#> 2  153.0046  -4.7089  POINT (5636527 -19294974)
#> 3  104.3261  -6.7541                POINT EMPTY
#> 4  124.9019   4.7817  POINT (11478868 18432292)
#> 5  126.7328   2.1643  POINT (11046583 19337712)
#> 6  153.2439  -5.6500  POINT (5589731 -19158700)
#> 7  142.8673  23.3882   POINT (6353660 16093116)
#> 8  152.6890  -5.5710  POINT (5673080 -19163103)

и вы можете написать и открыть с помощью QGIS: 

write_sf(sfdf_25833, 'mysf.gpkg')
...