0

I got a NetCDF file from the German Meteorological Service concerning mean temperatures in Europe (CDC FDP SERVER). The only thing I want to extract is the daily mean temperature for Bornholm, which is an island in the central Baltic.

I know how to extract information for certain coordinates (see code sample below). The only problem is that the file specific coordinates are 'rotated' which is why the geographic coordinates for Bornholm (extracted from GoogleMaps) are kind of useless.

packages <- c("RNetCDF",
              "ncdf4",
              "raster")

lapply(packages, require, character.only = TRUE)

x <- mean(14.68,15.16)        #coordinates for a rectangle around 
y <- mean(54.987,55.299)      #Bornholm extracted from GoogleMaps

temp <- nc_open("tas_decreg_europe_v20140120_20030101_20030131.nc")
temp

var <- ncvar_get(temp, "tas")
point <- var[x,y,]
as.data.frame(point)

To cut it short - Google uses a close variant of the Mercator projection. So how can I either convert the NetCDF file or the coordinates from GoogleMaps, so that I can find what I need. I could have bet that there's a simple solution out there but unfortunately not - at least I couldn't find one.

For information about the file generated by print(temp) see below:

File tas_decreg_europe_v20140120_20030101_20030131.nc (NC_FORMAT_CLASSIC):

     2 variables (excluding dimension variables):
        char rotated_pole[]   
            grid_mapping_name: rotated_latitude_longitude
            grid_north_pole_latitude: 39.25
            grid_north_pole_longitude: -162
        float tas[lon,lat,time]   
            long_name: Near-Surface Air Temperature
            units: K
            grid_mapping: rotated_pole
            _FillValue: 1.00000002004088e+20
            missing_value: 1.00000002004088e+20

     3 dimensions:
        lon  Size:1056
            standard_name: grid_longitude
            long_name: longitude
            units: degrees_east
            axis: X
        lat  Size:1026
            standard_name: grid_latitude
            long_name: latitude
            units: degrees_north
            axis: Y
        time  Size:31   *** is unlimited ***
            standard_name: time
            units: days since 2003-01-01 00:00:00
            calendar: standard

Any help is appreciated. Thanks a lot...

Robert
  • 133
  • 10

1 Answers1

0

You load the raster package, but you do not use it. Have you tried something like the below?

library(raster)  
x <- mean(14.68,15.16)  
y <- mean(54.987,55.299)
temp <- brick("tas_decreg_europe_v20140120_20030101_20030131.nc", var='tas')
extract(temp, cbind(x,y))
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Definitely the correct way to extract the data, but still shows NAs only. No wonder, either, if the extent cannot contain my coordinates (`extent : -28.342, 18.122, -23.342, 21.802 (xmin, xmax, ymin, ymax)`). So, seems that I still need to convert something as mentioned already. Any further recommendations...? – Robert Aug 24 '16 at 20:16
  • These coordinates are clearly not the standard coordinates for anywhere in Europe, so indeed they seem to be using a different origin than (Greenwich, Equator), in other words, a rotation (or perhaps an projection?). But I cannot guess what that rotation is if you do not tell us. There is probably some documentation on their website. Or perhaps there is info when you do `print(temp)` – Robert Hijmans Aug 25 '16 at 08:40
  • That's why I provided the LINK to the file. Anyway, I just added the information generated by print(temp) to the original question (see above). The data set description from the website says pretty much the same. Cheers – Robert Aug 25 '16 at 10:01
  • You did, but I cannot download it from where I am (nor should I have to). – Robert Hijmans Aug 25 '16 at 12:39