Do Era-interim datasets have wrong proj4strings?

Question

I'm working with Era-interim datasets. I'd like to extract weather data for some cities. The code and data are updated to github.

Firstly, I use raster to read in the file I downloaded from the website:

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

Then I load my cities:

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

... and convert them to an sf object:

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)

As windspeed and capital_sf have different CRS, I need to perform some transformation:

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)

Strangely, the proj4string changes, but the coordinates don't change.

To see if my transformation is successful, I make a plot:

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

here is the plot:

The longitude spans from -0.375 to 359.627 rather than from -180 to 180. Therefore, all cities in the eastern hemisphere are correctly marked, but all cities in the western hemisphere are missing.

I'm confused. Why doesn't st_transform work? Do I pass a wrong proj4string, or the function simply can't process customized CRS?

It's hard to help without a [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) but in the mean time, you could try sp::spTransform as an alternative to pinpoint your issue. — Michael Bird, Aug 10 '18 at 12:39
Thank you very much. I tried sp::spTransform. It produced the same results. I uploaded the code and data to https://github.com/yutuotuo84/era-interim-example It should work as a reproducible example. — fishinweb, Aug 10 '18 at 15:25

score 1 · Answer 1 · answered Aug 10 '18 at 20:59

This is a good reference concerning the format of the ERA-Interim dataset:

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

It specifically reports that:

Longitudes range from 0 to 360, which is equivalent to -180 to +180 in Geographic coordinate systems.

A quick and dirty way to get what you want could be to "move" the right side of your raster to the left side, then manually adjusting the extent so that it spans -180 to 180.
In that way, the raster is in a "standard GCS" representation and you can easily use it for plotting.

For example:

# 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)

, which seems more or less correct.

HTH!

Actually I want to process a lot of GRIB files using the same capital dataset. So instead of swap the rasters, I swapped the capitals: capitals <- capitals %>% mutate(long = ifelse(long < -0.375, long + 360, long)). The map doesn't look good, but the extract function gives the same result as I swap the raster. — fishinweb, Aug 11 '18 at 04:02

Do Era-interim datasets have wrong proj4strings?

1 Answers1