【发布时间】：2019-01-18 01:00:48
【问题描述】：
我正在使用 Era-interim 数据集。我想提取一些城市的天气数据。代码和数据更新为github。
首先，我使用光栅读取从网站下载的文件：
library(raster)
windspeed <- raster("data/10m_wind_speed_19950101.grib")
windspeed
# class : RasterLayer
# dimensions : 241, 480, 115680 (nrow, ncol, ncell)
# resolution : 0.75, 0.75 (x, y)
# extent : -0.375, 359.625, -90.375, 90.375 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +a=6367470 +b=6367470 +no_defs
然后我加载我的城市：
load("capitals.RData")
head(capitals)
# ID iso3 country capital long lat
# 1 1 AUS Australia Canberra 149.13 -35.31
# 2 2 AUT Austria Vienna 16.37 48.22
# 3 3 BEL Belgium Brussels 4.33 50.83
# 4 4 BGR Bulgaria Sofia 23.31 42.69
# 5 5 BRA Brazil Brasilia -47.91 -15.78
# 6 6 CAN Canada Ottawa -75.71 45.42
...并将它们转换为 sf 对象：
library(sf)
capitals_sf <- st_as_sf(capitals, coords = c("long", "lat"), crs = 4326)
capitals_sf
# Simple feature collection with 40 features and 4 fields
# geometry type: POINT
# dimension: XY
# bbox: xmin: -99.14 ymin: -35.31 xmax: 149.13 ymax: 60.17
# epsg (SRID): 4326
# proj4string: +proj=longlat +datum=WGS84 +no_defs
# First 10 features:
# ID iso3 country capital geometry
# 1 1 AUS Australia Canberra POINT (149.13 -35.31)
# 2 2 AUT Austria Vienna POINT (16.37 48.22)
# 3 3 BEL Belgium Brussels POINT (4.33 50.83)
# 4 4 BGR Bulgaria Sofia POINT (23.31 42.69)
# 5 5 BRA Brazil Brasilia POINT (-47.91 -15.78)
# 6 6 CAN Canada Ottawa POINT (-75.71 45.42)
# 7 7 CHN China Beijing POINT (116.4 39.93)
# 9 8 CYP Cyprus Nicosia POINT (33.38 35.16)
# 11 9 CZE Czech Republic Prague POINT (14.43 50.08)
# 12 10 DEU Germany Berlin POINT (13.38 52.52)
由于windspeed 和capital_sf 有不同的CRS，我需要进行一些转换：
newcrs <- crs(windspeed, asText=TRUE)
capitals_tf <- st_transform(capitals_sf, newcrs)
capital_tf
# Simple feature collection with 40 features and 4 fields
# geometry type: POINT
# dimension: XY
# bbox: xmin: -99.14 ymin: -35.31 xmax: 149.13 ymax: 60.17
# epsg (SRID): NA
# proj4string: +proj=longlat +a=6367470 +b=6367470 +no_defs
# First 10 features:
# ID iso3 country capital geometry
# 1 1 AUS Australia Canberra POINT (149.13 -35.31)
# 2 2 AUT Austria Vienna POINT (16.37 48.22)
# 3 3 BEL Belgium Brussels POINT (4.33 50.83)
# 4 4 BGR Bulgaria Sofia POINT (23.31 42.69)
# 5 5 BRA Brazil Brasilia POINT (-47.91 -15.78)
# 6 6 CAN Canada Ottawa POINT (-75.71 45.42)
# 7 7 CHN China Beijing POINT (116.4 39.93)
# 9 8 CYP Cyprus Nicosia POINT (33.38 35.16)
# 11 9 CZE Czech Republic Prague POINT (14.43 50.08)
# 12 10 DEU Germany Berlin POINT (13.38 52.52)
奇怪的是，proj4string 改变了，但是坐标没有改变。
为了看看我的改造是否成功，我做了一个情节：
plot(windspeed)
plot(capitals_tf, col = "black", add = TRUE)
剧情如下：
经度范围从 -0.375 到 359.627，而不是从 -180 到 180。因此，东半球的所有城市都被正确标记，但西半球的所有城市都缺失。
我很困惑。为什么st_transform 不起作用？我是否传递了错误的 proj4string，或者该函数根本无法处理自定义的 CRS？
【问题讨论】：
-
如果没有 reproducible example 就很难提供帮助，但与此同时，您可以尝试使用 sp::spTransform 作为替代方法来查明您的问题。
-
非常感谢。我试过 sp::spTransform。它产生了相同的结果。我将代码和数据上传到github.com/yutuotuo84/era-interim-example 它应该可以作为可重现的示例。
标签： r geospatial raster sf