R Shapefile未正确绘制纬度/经度

huangapple go评论113阅读模式
英文:

R Shapefile not plotting lat/lon correctly

问题

我在绘制一个形状文件时遇到了问题,似乎已经识别了几何和投影,并且我可以在geometry$num中看到适当的值(纬度约为-34,经度约为153),但在绘制时完全错误。我认为投影在某个地方出错了。

可以在此处访问形状文件:https://github.com/HaydenSchilling/Example_code_and_data/tree/master/shapefile

当前用于重现问题的代码是:

  1. library(sf)
  2. library(tidyverse)
  3. dat <- sf::read_sf("227151 (1)99%fulltrack_WGS84.shp")
  4. ggplot(dat) + geom_sf()

目前的绘图结果如下:
R Shapefile未正确绘制纬度/经度

如何手动修复投影?

英文:

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:

  1. library(sf)
  2. library(tidyverse)
  3. dat &lt;- sf::read_sf(&quot;227151 (1)99%fulltrack_WGS84.shp&quot;)
  4. ggplot(dat) + geom_sf()

The result plot currently looks like:
R Shapefile未正确绘制纬度/经度

How do I manually fix the projection?

答案1

得分: 4

以下是您要求的代码部分的翻译:

  1. # transform to original CRS (GDA94 / Geoscience Australia Lambert)
  2. dat_gda94 <- st_transform(dat_wgs84, dat_crs)
  3. dat_gda94
  4. #> Simple feature collection with 1 feature and 2 fields
  5. #> Geometry type: POLYGON
  6. #> Dimension: XY
  7. #> Bounding box: xmin: 1456105 ymin: -4173843 xmax: 1715909 ymax: -3755963
  8. #> Projected CRS: GDA94 / Geoscience Australia Lambert
  9. #> # A tibble: 1 × 3
  10. #> timestamp sn geometry
  11. #> * <chr> <chr> <POLYGON [m]>
  12. #> 1 2022-03-25-2022-06-18 227151 (1) ((1486283 -3925744, 1486170 -3926615, 148575…
  13. ggplot(dat_gda94) + geom_sf()

希望这对您有所帮助。

英文:

> (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:

  1. library(sf)
  2. #&gt; Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
  3. library(ggplot2)
  4. shp_url &lt;- &quot;https://github.com/HaydenSchilling/Example_code_and_data/raw/master/shapefile/227151%20(1)99%25fulltrack_WGS84.shp&quot;
  5. dat &lt;- sf::read_sf(paste0(&quot;/vsicurl/&quot;, shp_url))
  6. # let&#39;s check the object, coordinates and proj4str
  7. # X - easting, longitude
  8. # Y - northing, latitude
  9. dat
  10. #&gt; Simple feature collection with 1 feature and 2 fields
  11. #&gt; Geometry type: POLYGON
  12. #&gt; Dimension: XY
  13. #&gt; Bounding box: xmin: -36.21875 ymin: 150.2143 xmax: -32.31378 ymax: 152.6772
  14. #&gt; Projected CRS: GDA94 / Geoscience Australia Lambert
  15. #&gt; # A tibble: 1 &#215; 3
  16. #&gt; timestamp sn geometry
  17. #&gt; &lt;chr&gt; &lt;chr&gt; &lt;POLYGON [m]&gt;
  18. #&gt; 1 2022-03-25-2022-06-18 227151 (1) ((-34.00546 150.2143, -34.01342 150.2143, -3…
  19. st_coordinates(dat)[1:5,]
  20. #&gt; X Y L1 L2
  21. #&gt; [1,] -34.00546 150.2143 1 1
  22. #&gt; [2,] -34.01342 150.2143 1 1
  23. #&gt; [3,] -34.04276 150.2143 1 1
  24. #&gt; [4,] -34.04364 150.2143 1 1
  25. #&gt; [5,] -34.06255 150.2143 1 1
  26. dat_crs &lt;- st_crs(dat)
  27. dat_crs$proj4string
  28. #&gt; [1] &quot;+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&quot;
  29. # does not look like projected easting/northing ...
  30. # let&#39;s make an assumption based on a file name and set (NOT transform) it to WGS84 :
  31. dat_wgs84 &lt;- st_set_crs(dat, &quot;WGS84&quot;)
  32. #&gt; Warning: st_crs&lt;- : replacing crs does not reproject data; use st_transform for
  33. #&gt; that
  34. # swap lat/lon in geometry, g[,c(1,2)] &lt;- g[,c(2,1)] :
  35. dat_wgs84$geometry[[1]][[1]][,c(1,2)] &lt;- dat_wgs84$geometry[[1]][[1]][,c(2,1)]
  36. # and update bbox :
  37. attr(st_geometry(dat_wgs84), &quot;bbox&quot;) &lt;- sfheaders::sf_bbox(dat_wgs84)
  38. # check the result
  39. dat_wgs84
  40. #&gt; Simple feature collection with 1 feature and 2 fields
  41. #&gt; Geometry type: POLYGON
  42. #&gt; Dimension: XY
  43. #&gt; Bounding box: xmin: 150.2143 ymin: -36.21875 xmax: 152.6772 ymax: -32.31378
  44. #&gt; Geodetic CRS: WGS 84
  45. #&gt; # A tibble: 1 &#215; 3
  46. #&gt; timestamp sn geometry
  47. #&gt; * &lt;chr&gt; &lt;chr&gt; &lt;POLYGON [&#176;]&gt;
  48. #&gt; 1 2022-03-25-2022-06-18 227151 (1) ((150.2143 -34.00546, 150.2143 -34.01342, 15…
  49. mapview::mapview(dat_wgs84)

R Shapefile未正确绘制纬度/经度<!-- -->

  1. # transform to original CRS (GDA94 / Geoscience Australia Lambert)
  2. dat_gda94 &lt;- st_transform(dat_wgs84, dat_crs)
  3. dat_gda94
  4. #&gt; Simple feature collection with 1 feature and 2 fields
  5. #&gt; Geometry type: POLYGON
  6. #&gt; Dimension: XY
  7. #&gt; Bounding box: xmin: 1456105 ymin: -4173843 xmax: 1715909 ymax: -3755963
  8. #&gt; Projected CRS: GDA94 / Geoscience Australia Lambert
  9. #&gt; # A tibble: 1 &#215; 3
  10. #&gt; timestamp sn geometry
  11. #&gt; * &lt;chr&gt; &lt;chr&gt; &lt;POLYGON [m]&gt;
  12. #&gt; 1 2022-03-25-2022-06-18 227151 (1) ((1486283 -3925744, 1486170 -3926615, 148575…
  13. ggplot(dat_gda94) + geom_sf()

R Shapefile未正确绘制纬度/经度<!-- -->

<sup>Created on 2023-08-09 with reprex v2.0.2</sup>

huangapple
  • 本文由 发表于 2023年8月9日 10:04:28
  • 转载请务必保留本文链接:https://go.coder-hub.com/76864137-2.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定