Era-interim 数据集是不是有错误的 proj4strings？

Posted 2023-03-23

【中文标题】Era-interim 数据集是不是有错误的 proj4strings？

【英文标题】：Do Era-interim datasets have wrong proj4strings?Era-interim 数据集是否有错误的 proj4strings？

【发布时间】：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 它应该可以作为可重现的示例。

【参考方案1】：

这是关于 ERA-Interim 数据集格式的一个很好的参考：

https://confluence.ecmwf.int/display/CKB/ERA-Interim%3A+What+is+the+spatial+reference

它特别报告：

经度范围从 0 到 360，相当于地理坐标系中的 -180 到 +180。

获得所需内容的一种快速而肮脏的方法可能是“移动”栅格的右侧 到左侧，然后手动调整范围，使其跨越 -180 到 180。 这样，栅格就处于“标准 GCS”表示中，您可以轻松地 用它来绘图。

例如：

# create temporary raster, then "move" the right side to the left
tmp <- windspeed
tmp[, 1:240] <- windspeed[, 241:480]
tmp[, 241:480] <- windspeed[, 1:240]

# put data back in windspeed (not really needed) and update extent
windspeed <- tmp
extent(windspeed)@xmin <- extent(windspeed)@xmin -180
extent(windspeed)@xmax <- extent(windspeed)@xmax -180

windspeed
class       : RasterLayer 
dimensions  : 241, 480, 115680  (nrow, ncol, ncell)
resolution  : 0.75, 0.75  (x, y)
extent      : -180.375, 179.625, -90.375, 90.375  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +a=6367470 +b=6367470 +no_defs 
data source : in memory
names       : X10m_wind_speed_19950101 
values      : 0.9062432, 14.906  (min, max)

# now plot: 
capitals_sf <- st_as_sf(capitals, coords = c("long", "lat"), crs = 4326)

plot(windspeed)
plot(capitals_sf, col = "black", add = TRUE)

，这似乎或多或少是正确的。

HTH！

【讨论】：

非常感谢！你的回答启发了我！

其实我想用同一个大写数据集来处理很多GRIB文件。因此，我没有交换栅格，而是交换了大写字母：大写字母 % mutate(long = ifelse(long