1

I have yearly gridded temperature observation for all countries, and I picked up recent 35 years' data for climate research. However, in original data, yearly temperature observation has been stored by global level, and I intend to retrieve all coordinates with its metadata (here its metadata is 12 month' average monthly temperature) only for Germany. I used data.table package to read those data in a clean format in R. Now I need to retrieve all points with its metadata that falls in German polygon. I don't need to visualize them first because I have to compute annual temperature for each coordinate that falls in Germany. I used maptools, sp packages to get this done, but my output was not what I expected.

Here is how I load my data in R (old post can be found here: old post):

getwd()
setwd("stella/data/air_temp_1980_2014/")

tempdat <- list.files(path = getwd(), recursive = TRUE)
TempDatas <- lapply(tempdat, function(x) {
    data.table::fread(x, sep = " ", header = FALSE,
                      col.names = c("Long", "Lat", "Jan","Feb", "Mar", "April", "May",
                                    "Jun", "Jul", "Aug", "Sept", "Oct", "Nov", "Dec"))
})

Here is how data look like (below example data is yearly temperature observation for the whole world):

> head(temperaturesDat)
$air_temp.1980
          Long    Lat   Jan   Feb   Mar April   May   Jun   Jul   Aug  Sept   Oct   Nov   Dec
    1: -179.75  71.25 -24.7 -24.0 -23.2 -18.3  -8.4   0.2   1.5   0.6  -2.0  -9.6 -18.6 -22.4
    2: -179.75  68.75 -27.0 -28.2 -27.2 -21.6  -9.0   0.6   2.8   1.9  -0.2 -11.9 -22.7 -25.1
    3: -179.75  68.25 -27.8 -28.5 -27.5 -22.0  -9.5   0.4   3.0   1.8  -0.8 -12.7 -23.6 -26.8
    4: -179.75  67.75 -26.8 -26.6 -25.7 -20.5  -8.0   2.7   6.0   4.0   0.5 -12.2 -23.2 -27.3
    5: -179.75  67.25 -29.1 -28.4 -27.5 -22.3  -9.7   2.2   6.2   3.3  -1.3 -15.4 -26.4 -31.1
   ---                                                                                       
85790:  179.75 -87.75 -26.7 -39.0 -52.0 -54.9 -54.3 -55.8 -58.0 -59.8 -58.0 -49.6 -32.9 -24.9
85791:  179.75 -88.25 -27.4 -39.8 -53.1 -56.0 -55.5 -56.9 -59.1 -61.2 -59.2 -50.7 -33.4 -25.3
85792:  179.75 -88.75 -27.5 -40.3 -53.7 -56.6 -56.1 -57.5 -59.8 -61.9 -59.8 -51.2 -33.5 -25.4
85793:  179.75 -89.25 -27.4 -40.3 -53.9 -56.7 -56.1 -57.7 -59.9 -62.2 -59.8 -51.3 -33.4 -25.1
85794:  179.75 -89.75 -27.2 -40.2 -54.0 -56.4 -55.8 -57.8 -59.7 -62.4 -59.8 -51.2 -33.5 -24.9

$air_temp.1981
          Long    Lat   Jan   Feb   Mar April   May   Jun   Jul   Aug  Sept   Oct   Nov   Dec
    1: -179.75  71.25 -22.4 -23.7 -22.2 -16.4  -4.2   0.7   1.7   0.7  -3.1  -9.0 -16.2 -22.0
    2: -179.75  68.75 -26.6 -26.5 -25.5 -16.8  -4.2   0.7   1.8   1.2  -2.5  -8.5 -17.6 -26.9
    3: -179.75  68.25 -27.1 -27.3 -26.1 -17.2  -5.0   0.7   2.1   1.1  -3.1  -9.8 -18.7 -27.3
    4: -179.75  67.75 -25.7 -26.2 -24.8 -15.5  -3.7   3.0   5.3   3.5  -1.6  -9.7 -18.4 -26.1
    5: -179.75  67.25 -27.5 -28.6 -27.3 -17.2  -5.5   2.5   5.7   3.0  -3.3 -13.2 -21.7 -28.6
   ---                                                                                       
85790:  179.75 -87.75 -25.7 -38.3 -50.0 -57.7 -53.3 -55.6 -54.7 -54.3 -60.2 -53.5 -39.6 -27.4
85791:  179.75 -88.25 -25.8 -39.0 -51.3 -58.7 -54.2 -56.8 -55.7 -55.1 -61.3 -54.7 -40.5 -28.0
85792:  179.75 -88.75 -25.6 -39.3 -51.9 -59.1 -54.5 -57.5 -56.3 -55.4 -61.9 -55.4 -40.9 -28.3
85793:  179.75 -89.25 -25.2 -39.3 -52.2 -59.1 -54.2 -57.8 -56.3 -55.3 -62.0 -55.6 -41.1 -28.1
85794:  179.75 -89.75 -24.8 -39.0 -52.3 -58.6 -53.8 -57.9 -56.0 -55.1 -62.0 -55.6 -41.4 -28.1

Here is what I tried to make this happen:

library(rgeos)
library(maptools)
library(sp)

data(wrld_simpl)
germany<- wrld_simpl[wrld_simpl@data$NAME =="Germany",]
pointDatas <- lapply(TempDatas, function(x) {
  pts <- data.frame(Long=x$Long, Lat=x$Lat)
  res <- pts[, c(1,2)]
})

then I tried Map and SpatialPointsDataFrame to get the coordinates that fall in Germany:

co<-proj4string(germany)
points_in_germany <- Map(SpatialPointsDataFrame, coords=pointDatas, 
                         data=pointDatas, proj4string = CRS(co))

I tried my approach but I got an error down below:

Error in dots[[3L]][[1L]] : this S4 class is not subsettable
In addition: There were 15 warnings (use warnings() to see them)

but the result does not seem what I want. How can I improve implementation above? Any way to make this happen in R? Any idea?

Desired output:

Basically, I need to retrieve all points with its metadata (12-month' temperature observations) that fall in Germany's polygon in new data.frame, then compute annual temperature for each point and add it as a new column. How can I get this done? Any way to do this more efficiently in R? Any thought?

Update:

How to merge multiple SpatialPointsDataFrame into one? Any idea?

  • 1
    I went ahead and addressed your question but, you did not at all follow my advice from your previous post. I explained why the structure of your data required a workaround, which I even provided code for. If you are going to ignore our recommendations, then why ask? – Jeffrey Evans Mar 30 '18 at 21:08
  • It is loading just fine for me with read.table. What is not working? You could just coerce the data.table to a data.frame using as.data.frame(). – Jeffrey Evans Mar 30 '18 at 21:21
  • @JeffreyEvans I tried to subset and merged to SpatialPointsDataFrame, but I got an error. I tried like this: merge(germany.pts, stations[stations$year >= 1980 & stations$year <=2014,], by = "ID"), because I need continuous 35 years' temperature data. Plus, I tried also this: merge(germany.pts, stations[data.table::between(stations$year, 1980, 2014)], by = "ID"). How can I get correct time ranges for merging ? Thank you very much. – Dan Mar 30 '18 at 22:28
  • This is because you have 181 stations in Germany and the data is ordered so, there should be 181 observations per year so, the dimensions will not match. You need to put some though in how you want to analyze this data because having a duplicate spatial feature, for a fixed location, for every year is not good practice nor very useful. It would be best to preform analysis on the subset data and then merge results back to the station locations. – Jeffrey Evans Mar 30 '18 at 22:51
4

Following your question on reading and formatting this specific climate data:

First, add libraries and country boundaries.

library(sp)
library(raster)
library(maptools) 

data(wrld_simpl)
germany<- wrld_simpl[wrld_simpl@data$NAME =="Germany",]

Now, assuming we have the data.frame "stations" created in your previous question (at this point, not sure why you split it into a list) we then create a unique station identifier using a concatenated vector of the coordinate pairs that will be used to match back to the climate data. Then we create a data.frame containing these unique ids along with coordinate pairs [x,y]. This data.frame is then coerced into a SpatialPointsDataFrame object.

xy <- unique( paste(stations[,1], stations[,2], sep="_")  )
stations.xy <- data.frame(long = as.numeric(unlist(lapply(strsplit(xy, "_"), function(x) {x[1]}))), 
                          lat =  as.numeric(unlist(lapply(strsplit(xy, "_"), function(x) {x[2]}))))
    stations.xy$ID <- xy                      
    coordinates(stations.xy) <- ~long+lat

For sub-setting to the Germany polygon boundary, first we must match the proj4string (assuming that the data are in the same projections). Then we can use a spatial overlay function to subset the station locations.

proj4string(stations.xy) <- proj4string(germany)
germany.pts <- raster::intersect(stations.xy, germany)

plot(germany)
  plot(germany.pts,pch=19,add=TRUE)

We can then match the concatenated ids to create an bracket index query that will subset the full "stations" dataset to observations only occurring within the polygon boundary of Germany.

idx <- which(stations$ID %in% germany.pts$ID)
  stations <- stations[idx,]

Now, you have a data.frame that matches your station locations. It can easily be subset and merged to the germany.pts SpatialPointsDataFrame.

stations.1929 <- merge(germany.pts, stations[stations$year == 1929,], by = "ID")
  head(stations.1929@data)
  • 1
    If you are going to change code then, you are on your own. With the data you provided, the solution in the previous post works just fine. Not exactly sure what is not "desirable" about the result in the first part of the code. It is retaining data integrity, returning ordered monthlies and providing a data structure to allow for coercion to an sp class points object with linkage back to the full yearly climate data. – Jeffrey Evans Mar 30 '18 at 21:19
  • thank you very much for your help, it worked well except a bit slow when I used bigger data. Perhaps not using for loop. How can I refine your solution even that work fast for big climate dataset? Any hint? – Dan Mar 30 '18 at 21:33
  • 1
    Here is where you can start playing with "apply" type functions along with multi-treading. Since you have a general framework you can also start replacing some the these functions with functions from the tidyverse packages. – Jeffrey Evans Mar 30 '18 at 21:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.