投影到 R 中的 EPSG:25833 时 QGIS 中的未知 CRS
Unknown CRS in QGIS when projecting to EPSG:25833 in R
我想在 R 中将空间数据帧投影到 EPSG 25833 ,但 QGIS 似乎不知道它(为了重现性,我使用在 [=24 中创建的代码 jazzurro =] 对 的回答稍有改动)
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
现在,当我将 shapefile 加载到 QGIS 时,它不知道投影。
有什么想法吗?
让你的 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"))`
您是在告诉数据框您的点所在的 crs,因此您应该指定 4326,因为您的原始数据是 lon/lat。
所以应该是:
spdf <- SpatialPointsDataFrame(coords = xy, data = mydf,
proj4string = CRS("+proj=longlat +datum=WGS84"))
然后您可以使用 spTransform
:
将数据转换为另一个 CRS
spTransform(spdf, CRS('+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'))
对于此特定数据,我们收到错误消息,因为其中一个点未转换为您的目标 CRS:
Error in spTransform(xSP, CRSobj, ...) : failure in points 3 In
addition: Warning message: In spTransform(xSP, CRSobj, ...) : 1
projected point(s) not finite
我更喜欢在 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')
我想在 R 中将空间数据帧投影到 EPSG 25833 ,但 QGIS 似乎不知道它(为了重现性,我使用在 [=24 中创建的代码 jazzurro =] 对
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
现在,当我将 shapefile 加载到 QGIS 时,它不知道投影。
有什么想法吗?
让你的 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"))`
您是在告诉数据框您的点所在的 crs,因此您应该指定 4326,因为您的原始数据是 lon/lat。
所以应该是:
spdf <- SpatialPointsDataFrame(coords = xy, data = mydf,
proj4string = CRS("+proj=longlat +datum=WGS84"))
然后您可以使用 spTransform
:
spTransform(spdf, CRS('+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'))
对于此特定数据,我们收到错误消息,因为其中一个点未转换为您的目标 CRS:
Error in spTransform(xSP, CRSobj, ...) : failure in points 3 In addition: Warning message: In spTransform(xSP, CRSobj, ...) : 1 projected point(s) not finite
我更喜欢在 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')