如何将SpatRaster的比例尺单位从km更改为mi?

s4n0splo  于 2023-02-01  发布在  其他
关注(0)|答案(1)|浏览(194)

我正在尝试使用SpatRaster底图在ggplot中创建Map(CRS 26913-NAD83,UTM ZONE 13N)。当我尝试使用函数annotation_scale添加比例尺时,scale_bar所表示的距离以千米为单位。但是,当我尝试将plot_unit参数更改为"mi"表示英里时,返回的Map仍有一个以公里为单位的不准确比例尺。我正在尝试获得以英里为单位的准确比例尺。
我有一个名为img的NAIP影像的SpatRaster,我正尝试使用ggplot绘制它,以便添加其他sf对象。我无法将比例尺从km更改为mi。我尝试了三种不同的函数,包括annotation_scaleggsn::scalebarterra::sbar。它当前位于坐标参考系NAD 83 UTM ZONE 13N(EPSG:第26913号)。

> class(img) 
[1] "SpatRaster"attr(,"package")
[1] "terra"
> crs(img)
[1] "PROJCRS["NAD83 / UTM zone 13N",\n    BASEGEOGCRS["NAD83",\n        DATUM["North American Datum 1983",\n            ELLIPSOID["GRS 1980",6378137,298.257222101,\n                LENGTHUNIT["metre",1]]],\n        PRIMEM["Greenwich",0,\n            ANGLEUNIT["degree",0.0174532925199433]],\n        ID["EPSG",4269]],\n    CONVERSION["UTM zone 13N",\n        METHOD["Transverse Mercator",\n            ID["EPSG",9807]],\n        PARAMETER["Latitude of natural origin",0,\n            ANGLEUNIT["degree",0.0174532925199433],\n            ID["EPSG",8801]],\n        PARAMETER["Longitude of natural origin",-105,\n            ANGLEUNIT["degree",0.0174532925199433],\n            ID["EPSG",8802]],\n        PARAMETER["Scale factor at natural origin",0.9996,\n            SCALEUNIT["unity",1],\n            ID["EPSG",8805]],\n        PARAMETER["False easting",500000,\n            LENGTHUNIT["metre",1],\n            ID["EPSG",8806]],\n        PARAMETER["False northing",0,\n            LENGTHUNIT["metre",1],\n            ID["EPSG",8807]]],\n    CS[Cartesian,2],\n        AXIS["(E)",east,\n            ORDER[1],\n            LENGTHUNIT["metre",1]],\n        AXIS["(N)",north,\n            ORDER[2],\n            LENGTHUNIT["metre",1]],\n    USAGE[\n        SCOPE["Engineering survey, topographic mapping."],\n        AREA["North America - between 108°W and 102°W - onshore and offshore. Canada - Northwest Territories; Nunavut; Saskatchewan. United States (USA) - Colorado; Montana; Nebraska; New Mexico; North Dakota; Oklahoma; South Dakota; Texas; Wyoming."],\n        BBOX[28.98,-108,84,-102]],\n    ID["EPSG",26913]]"

在这里,我尝试使用annotation_scaleMapSpatRaster。

# Map with annotation_scale with default settings (scalebar in km)
sat_1 <- ggplot() +
    geom_spatraster_rgb(data = img) +
    annotation_scale(location = "bl") +
    ggtitle("sat_1")

sat_1

Map with annotation_scale and default arguments (scale bar is accurate)

# Map with annotation_scale with argument "mi"
sat_2 <- ggplot() +
    geom_spatraster_rgb(data = img) +
    annotation_scale(location = "bl", plot_unit = "mi") +
    ggtitle("sat_2")

sat_2

Map with annotation_scale with plot_unit argument of miles (scalebar is wildly off)

# Map with ggsn::scalebar, didn't even work - see Error message below:
sat_3 <- ggplot() +
    geom_spatraster_rgb(data = img) +
    ggsn::scalebar(data = img, dist = 1, location = "bottomleft", transform = FALSE) +
    ggtitle("sat_3")

sat_3

SpatRaster resampled to ncells = 500379
Error: [subset] invalid name(s)

关于如何使用UTM或其他投影获得以英里为单位的精确比例尺,有人有什么建议吗?

qaxu7uf2

qaxu7uf21#

你可以用ggsn来做,问题是ggsn::scalebar()需要一个与你的SpatRaster相同范围的data对象。
您可以从SpatRaster的范围创建一个SpatVector(确保两者具有相同的CRS),然后使用sf::st_as_sf()将其传递给sf。
现在只需在方便的时候使用dist_unit参数:

library(tidyterra)
library(ggplot2)
library(terra)

# Get the data from extent using maptiles
x <- rast(xmin = -103.50, xmax = -103.38, ymin = 29.58, ymax = 29.70) %>%
  project("EPSG: 26913")

img <- maptiles::get_tiles(x, "Esri.WorldImagery") %>%
  project(pull_crs(x)) %>%
  crop(x)

# This is for checking with terra
# For locating
ext <- as.vector(ext(img))[c(1, 3)]

plotRGB(img, mar = c(3, 0, 0, 0))
sbar(6000, type = "bar", xy = ext - c(0, 1200), label = c("0", "3", "6 km"))

# On miles
sbar(6000 / 0.621371, type = "bar", xy = ext - c(0, 1200 * 2), label = c("0", "3", "6 mi"))

# (Sorry, image seems to be cropped with hop you get the idea)

# With ggsn we need a sf object
v <- ext(img) %>% as.polygons()
crs(v) <- pull_crs(img)

ggplot() +
  geom_spatraster_rgb(data = img) +
  ggsn::scalebar(
    data = sf::st_as_sf(v),
    location = "bottomleft",
    dist = 3,
    st.size = 4,
    st.dist = 0.05,
    transform = FALSE,
    dist_unit = "km"
  ) +
  coord_sf() +
  labs(
    x = "",
    y = ""
  )

# On miles
ggplot() +
  geom_spatraster_rgb(data = img) +
  ggsn::scalebar(
    data = sf::st_as_sf(v),
    location = "bottomleft",
    dist = 3,
    st.size = 4,
    st.dist = 0.05,
    transform = FALSE,
    dist_unit = "mi"
  ) +
  coord_sf() +
  labs(
    x = "",
    y = ""
  )

创建于2023年1月28日,使用reprex v2.0.2

相关问题