英文:
R Shapefile not plotting lat/lon correctly
问题
我在绘制一个shapefile时遇到了问题,似乎已经识别出了几何和投影,并且在geometry$num中可以看到适当的值(纬度约为-34,经度约为153),但是绘制出来的结果完全错误。我猜测投影在某个地方出错了。
可以在这里访问shapefile:https://github.com/HaydenSchilling/Example_code_and_data/tree/master/shapefile
目前用于重现问题的代码是:
library(sf)
library(tidyverse)
dat <- sf::read_sf("227151 (1)99%fulltrack_WGS84.shp")
ggplot(dat) + geom_sf()
我该如何手动修复投影?
英文:
I'm having trouble plotting a shapefile correctly, it seems to have the geometry and projection recognised and I can see appropriate values in geometry$num (lat ~-34, lon ~153) but when plotted it is completely wrong. I assume the projection is wrong somewhere.
Shapefile can be accessed here: https://github.com/HaydenSchilling/Example_code_and_data/tree/master/shapefile
current code to reproduce issue is:
library(sf)
library(tidyverse)
dat <- sf::read_sf("227151 (1)99%fulltrack_WGS84.shp")
ggplot(dat) + geom_sf()
The result plot currently looks like:
How do I manually fix the projection?
答案1
得分: 4
> (纬度 ~-34,经度~153)
看起来纬度和经度被交换了,而且没有进行投影坐标系的转换,只是替换了坐标参考系统。
让我们尝试修复这个问题:
library(sf)
#> 链接到 GEOS 3.9.3、GDAL 3.5.2、PROJ 8.2.1;sf_use_s2() 为 TRUE
library(ggplot2)
shp_url <- "https://github.com/HaydenSchilling/Example_code_and_data/raw/master/shapefile/227151%20(1)99%25fulltrack_WGS84.shp"
dat <- sf::read_sf(paste0("/vsicurl/", shp_url))
# 让我们检查对象、坐标和 proj4str
# X - 东经,经度
# Y - 北纬,纬度
dat
#> 简单要素集,有 1 个要素和 2 个字段
#> 几何类型: POLYGON
#> 维度: XY
#> 边界框: xmin: -36.21875 ymin: 150.2143 xmax: -32.31378 ymax: 152.6772
#> 投影坐标系: GDA94 / Geoscience Australia Lambert
#> # A tibble: 1 × 3
#> timestamp sn geometry
#> <chr> <chr> <POLYGON [m]>
#> 1 2022-03-25-2022-06-18 227151 (1) ((-34.00546 150.2143, -34.01342 150.2143, -3…
st_coordinates(dat)[1:5,]
#> X Y L1 L2
#> [1,] -34.00546 150.2143 1 1
#> [2,] -34.01342 150.2143 1 1
#> [3,] -34.04276 150.2143 1 1
#> [4,] -34.04364 150.2143 1 1
#> [5,] -34.06255 150.2143 1 1
dat_crs <- st_crs(dat)
dat_crs$proj4string
#> [1] "+proj=lcc +lat_0=0 +lon_0=134 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"
# 看起来不像是投影的东经/北纬...
# 让我们根据文件名做一个假设,并将其设置(而不是转换)为 WGS84:
dat_wgs84 <- st_set_crs(dat, "WGS84")
#> 警告: st_crs<- : replacing crs does not reproject data; use st_transform for
#> that
# 在几何图形中交换纬度/经度,g[,c(1,2)] <- g[,c(2,1)]:
dat_wgs84$geometry[[1]][[1]][,c(1,2)] <- dat_wgs84$geometry[[1]][[1]][,c(2,1)]
# 并更新边界框:
attr(st_geometry(dat_wgs84), "bbox") <- sfheaders::sf_bbox(dat_wgs84)
# 检查结果
dat_wgs84
#> 简单要素集,有 1 个要素和 2 个字段
#> 几何类型: POLYGON
#> 维度: XY
#> 边界框: xmin: 150.2143 ymin: -36.21875 xmax: 152.6772 ymax: -32.31378
#> 大地坐标系: WGS 84
#> # A tibble: 1 × 3
#> timestamp sn geometry
#> * <chr> <chr> <POLYGON [°]>
#> 1 2022-03-25-2022-06-18 227151 (1) ((150.2143 -34.00546, 150.2143 -34.01342, 15…
mapview::mapview(dat_wgs84)
<!-- -->
# 转换回原始的 CRS(GDA94 / Geoscience Australia Lambert)
dat_gda94 <- st_transform(dat_wgs84, dat_crs)
dat_gda94
#> 简单要素集,有 1 个要素和 2 个字段
#> 几何类型: POLYGON
#> 维度: XY
#> 边界框: xmin: 1456105 ymin: -4173843 xmax: 1715909 ymax: -3755963
#> 投影坐标系: GDA94 / Geoscience Australia Lambert
#> # A tibble: 1 × 3
#> timestamp sn geometry
#> * <chr> <chr> <POLYGON [m]>
#> 1 2022-03-25-2022-06-18 227151 (1) ((1486283 -3925744, 1486170 -3926615, 148575…
ggplot(dat_gda94) + geom_sf()
<!-- -->
<sup>使用 reprex v2.0.2 在 2023-08-09 创建</sup>
英文:
> (lat ~-34, lon ~153)
It looks like latitude & longitude have been swapped, and instead of transforming to projected CRS, CRS was just replaced.
Let's try to fix this:
library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(ggplot2)
shp_url <- "https://github.com/HaydenSchilling/Example_code_and_data/raw/master/shapefile/227151%20(1)99%25fulltrack_WGS84.shp"
dat <- sf::read_sf(paste0("/vsicurl/", shp_url))
# let's check the object, coordinates and proj4str
# X - easting, longitude
# Y - northing, latitude
dat
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: -36.21875 ymin: 150.2143 xmax: -32.31378 ymax: 152.6772
#> Projected CRS: GDA94 / Geoscience Australia Lambert
#> # A tibble: 1 × 3
#> timestamp sn geometry
#> <chr> <chr> <POLYGON [m]>
#> 1 2022-03-25-2022-06-18 227151 (1) ((-34.00546 150.2143, -34.01342 150.2143, -3…
st_coordinates(dat)[1:5,]
#> X Y L1 L2
#> [1,] -34.00546 150.2143 1 1
#> [2,] -34.01342 150.2143 1 1
#> [3,] -34.04276 150.2143 1 1
#> [4,] -34.04364 150.2143 1 1
#> [5,] -34.06255 150.2143 1 1
dat_crs <- st_crs(dat)
dat_crs$proj4string
#> [1] "+proj=lcc +lat_0=0 +lon_0=134 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"
# does not look like projected easting/northing ...
# let's make an assumption based on a file name and set (NOT transform) it to WGS84 :
dat_wgs84 <- st_set_crs(dat, "WGS84")
#> Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
#> that
# swap lat/lon in geometry, g[,c(1,2)] <- g[,c(2,1)] :
dat_wgs84$geometry[[1]][[1]][,c(1,2)] <- dat_wgs84$geometry[[1]][[1]][,c(2,1)]
# and update bbox :
attr(st_geometry(dat_wgs84), "bbox") <- sfheaders::sf_bbox(dat_wgs84)
# check the result
dat_wgs84
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 150.2143 ymin: -36.21875 xmax: 152.6772 ymax: -32.31378
#> Geodetic CRS: WGS 84
#> # A tibble: 1 × 3
#> timestamp sn geometry
#> * <chr> <chr> <POLYGON [°]>
#> 1 2022-03-25-2022-06-18 227151 (1) ((150.2143 -34.00546, 150.2143 -34.01342, 15…
mapview::mapview(dat_wgs84)
<!-- -->
# transform to original CRS (GDA94 / Geoscience Australia Lambert)
dat_gda94 <- st_transform(dat_wgs84, dat_crs)
dat_gda94
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 1456105 ymin: -4173843 xmax: 1715909 ymax: -3755963
#> Projected CRS: GDA94 / Geoscience Australia Lambert
#> # A tibble: 1 × 3
#> timestamp sn geometry
#> * <chr> <chr> <POLYGON [m]>
#> 1 2022-03-25-2022-06-18 227151 (1) ((1486283 -3925744, 1486170 -3926615, 148575…
ggplot(dat_gda94) + geom_sf()
<!-- -->
<sup>Created on 2023-08-09 with reprex v2.0.2</sup>
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论