geom_contour() 将 WGS84 坐标投影到 UTM 后出错

Question

我有测深数据（在 WGS84 Lat/Long 中），我想使用

ggplot

将其绘制为等高线。我需要 UTM 中的测深轮廓。使用 Lat/Long 坐标时，我可以获得轮廓，但使用投影 UTM 坐标时却不能。

require(marmap)
bathy <- getNOAA.bathy(lon1 = -156.4 , lon2 = -156.95, lat1 = 20.55, lat2 = 21.05 , resolution = 1)
bathy.df<-fortify(bathy)
bathy.df<-cbind(bathy.df,
 as.data.frame(terra::project(as.matrix(bathy.df[,c(1,2)]), "+proj=longlat", 
"+proj=utm +zone=4 +units=km")))
colnames(bathy.df)<-c("x", "y", "z", "x.utm", "y.utm")

require(ggplot2)
ggplot(data=bathy.df, aes(x=x, y=y, z=z)) +
  geom_contour(binwidth=100)

ggplot(data=bathy.df, aes(x=x.utm, y=y.utm, z=z))+geom_contour(binwidth=100)
Warning messages:
1: `stat_contour()`: Zero contours were generated 
2: In min(x) : no non-missing arguments to min; returning Inf
3: In max(x) : no non-missing arguments to max; returning -Inf

Answer 1

好吧，我明白了；发布答案以防其他人遇到同样的问题。实际上绘制 y~x 和 y.utm~x.utm 有助于可视化正在发生的事情。

ggplot(bathy.df, aes(x=x,y=y,color=z))+geom_point()

ggplot(bathy.df, aes(x=x.utm,y=y.utm,color=z))+geom_point()

基本上，它们并没有完美排列。在我将 UTM 坐标四舍五入后，我现在可以创建轮廓了。

Answer 2

发生的事情是

bathy

是

longlat

中的规则网格，但是当您将每个点（代表网格的中心）投影到不同的 crs (

x.utm, y.utm

) 时，这些点会发生变化，它们的当前位置也会发生变化不再符合规则网格，这是

geom_contour()

所必需的。 https://ggplot2.tidyverse.org/reference/geom_contour.html

这些函数需要常规数据，其中x和y坐标形成一个等距的网格，x和y的每个组合出现一次。

我的建议是不要使用

terra

来投影点，而是将

bathy

转换为 SpatRaster 网格，然后“投影”（或扭曲）SpatRaster，这是使用常规网格的正确方法。然后，作为建议（并且因为您正在使用

ggplot2

），您可以使用

tidyterra

来创建绘图轮廓。

转换只能发生在绘图阶段（使用

coord_sf()

）或创建一个新的投影 SpatRaster（见

bathy_terra_proj

）：

看一些例子：

library(marmap)
#> Registered S3 methods overwritten by 'adehabitatMA':
#>   method                       from
#>   print.SpatialPixelsDataFrame sp  
#>   print.SpatialPixels          sp
#> 
#> Attaching package: 'marmap'
#> The following object is masked from 'package:grDevices':
#> 
#>     as.raster
library(terra)
#> terra 1.7.23
library(tidyterra)
#> 
#> Attaching package: 'tidyterra'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(ggplot2)

bathy <- getNOAA.bathy(
  lon1 = -156.4, lon2 = -156.95,
  lat1 = 20.55, lat2 = 21.05,
  resolution = 1
)
#> Querying NOAA database ...
#> This may take seconds to minutes, depending on grid size
#> Building bathy matrix ...

# To terra

bathy_terra <- fortify(bathy) %>% rast(crs = "+proj=lonlat")

bathy_terra
#> class       : SpatRaster 
#> dimensions  : 30, 32, 1  (nrow, ncol, nlyr)
#> resolution  : 0.01612903, 0.01616379  (x, y)
#> extent      : -156.9081, -156.3919, 20.55754, 21.04246  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#> source(s)   : memory
#> name        :         z 
#> min value   : -478.0886 
#> max value   : 1519.8207

# Plot with ggplot2

gpl <- ggplot() +
  geom_spatraster_contour(
    data = bathy_terra,
    binwidth = 100
  )

gpl

# With coords utm
gpl +
  coord_sf(crs = "+proj=utm +zone=4 +units=km") +
  labs(title = "With utm",
       subtitle = "transformation happened in the plot, not in the raster")

# Variatios
## With contours and coords

ggplot() +
  geom_spatraster_contour_filled(data = bathy_terra, binwidth = 100) +
  coord_sf(crs = "+proj=utm +zone=4 +units=km")

# With labeled contours
# See https://dieghernan.github.io/tidyterra/articles/faqs.html#label-contour
# Need first to project since text contour needs a df

bathy_terra_proj <- bathy_terra %>% project("+proj=utm +zone=4 +units=km")

bathy_terra_proj
#> class       : SpatRaster 
#> dimensions  : 31, 31, 1  (nrow, ncol, nlyr)
#> resolution  : 1.732274, 1.732274  (x, y)
#> extent      : 717.3869, 771.0874, 2275.362, 2329.063  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=utm +zone=4 +datum=WGS84 +units=km +no_defs 
#> source(s)   : memory
#> name        :         z 
#> min value   : -455.5618 
#> max value   : 1285.8115

require("metR")
#> Loading required package: metR
library(metR)

ggplot(bathy_terra_proj) +
  geom_spatraster_contour(data = bathy_terra_proj, binwidth = 100) +
  geom_text_contour(aes(x, y, z = z),
                    binwidth = 100,
                    check_overlap = TRUE, stroke = 0.2,
                    stroke.colour = "white"
  ) +
  labs(x = "", y = "")

^{创建于 2023-04-14 与 reprex v2.0.2}