Я работаю с наборами данных 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?
Спасибо большое. Пробовал sp :: spTransform. Это дало те же результаты. Я загрузил код и данные в github.com/yutuotuo84/era-interim-example. Он должен работать как воспроизводимый пример.
Это хорошая справочная информация о формате набора данных ERA-Interim:
https://confluence.ecmwf.int/display/CKB/ERA-Interim%3A+What+is+the+spatial+reference
В нем конкретно говорится, что:
Longitudes range from 0 to 360, which is equivalent to -180 to +180 in Geographic coordinate systems.
Быстрый и грязный способ получить то, что вы хотите, - это "переместить" правую часть вашего растра. влево, затем вручную отрегулируйте экстент так, чтобы он составлял от -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, используя один и тот же набор данных с заглавными буквами. Поэтому вместо того, чтобы менять местами растры, я поменял местами заглавные буквы: capitals <- capitals%>% mutate (long = ifelse (long <-0.375, long + 360, long)). Карта выглядит не очень хорошо, но функция извлечения дает тот же результат, что и я меняю растр местами.
Трудно помочь без воспроизводимый пример, но пока вы можете попробовать sp :: spTransform в качестве альтернативы, чтобы точно определить вашу проблему.