domingo, 25 de noviembre de 2018

Plotting wind highways using rWind



Hi there!

Our manuscript about rWind R package has been recently accepted for publication in Ecography! As you know, rWind is a tool used to download and manage wind data, with some utilities that make easy to include wind information in ecological or evolutionary analyses (or others!).

Though there are several examples in the Supporting material of the publication, and others in the vignette included in the CRAN version, here you have a full tutorial with the main utilities of rWind... enjoy!

Computing wind connectivity of Azores Archipelago


In this example, we will compute wind connectivity of Azores Islands with different points of Europe, Africa and North America during June 2018. I will divide each task of this tutorial in different sections for clarification.

 

Downloading wind time series using lubridate and wind.dl_2()


In the last version of rWind (1.0.4), we can use lubridate to create a list of dates which will be used latter with wind.dl_2 to download wind time series. The object obtained with wind.dl_2 is an "rWind_series" "list". This object can be easily used with other rWind functions as we'll see later.

 # Make sure that you have the last version of rWind  
   
 devtools::install_github("jabiologo/rWind")  
 library(rWind)  
 library(lubridate)  
   
 # Here, we use ymd_hms from lubridate package to create a sequence of dates  
   
 dt <- seq(ymd_hms(paste(2018,6,1,12,00,00, sep="-")),  
      ymd_hms(paste(2018,6,30,12,00,00, sep="-")),by="1 days")  
   
 # Now we can use wind.dl_2 with this sequence of dates. Have into account  
 # that it could take a while, they are 30 datasets and it's a big area.   
   
 ww <- wind.dl_2(dt,-85,5,20,60)  


Obtaining maximum wind speed during the time series with dplyr


We can use the rWind function called "tidy" to play with the "rWind_series" "list" object to obtain, for example, the maximum speed reported for each cell in the study area during the time series.

 # After tidy our wind data, we can use dplyr utilities to easily  
 # obtain several metrics from our wind data. Notice that wind average  
 # is a special case, and it will be discussed later  
   
 t_ww <- tidy(ww)  
   
 library(dplyr)  
 g_ww <- t_ww %>% group_by(lat, lon)  
   
 # Now, we select only the maximum speed values for each coordinate  
   
 max_ww <- g_ww %>% summarise(speed = max(speed))  
 maxw <- cbind(max_ww$lon, max_ww$lat, max_ww$speed)  
 head(maxw)  
   
 # We can also check the maximum speed value in the whole study area   
   
 max(maxw[,3])  
   
 # Finally, we can transform this data into a raster layer  
   
 rmax <- rasterFromXYZ(maxw)  
 acol <- colorRampPalette(c("white", "blue", "darkblue"))  
 plot(rmax, col=acol(1000), main= "Maximum wind speed reported",  
     xlab="Longitude", ylab="Lattitude")  
   
 library(rworldmap)  
 lines(getMap(resolution = "high"), lwd=2)  



Computing averages of wind speeds and directions


Since wind data is a vectorial magnitude given in two components (U and V), to compute speed and direction averages we should use both vectors. To make it easier, rWind includes wind.mean function. Notice that wind.mean is specially developed to work with wind.dl_2 outputs, so it will require "rWind_series" object as input.
 # We can easily use wind.mean() function  
   
 w_mean <- wind.mean(ww)  
   
 # And then use wind2raster directly to create a raster layer  
   
 r_mean <- wind2raster(w_mean)  

 # We can plot a subset of this raster around Azores Islands.  
 # Using "arrowDir" we can include arrows for wind direction with   
 # "Arrowhead" function from "shape" R package  
   
 plot(r_mean$wind.speed, main= "Wind speed and direction average",
      col=acol(1000), xlim=c(-38,-18), ylim=c(33,43))  
 alpha <- arrowDir(w_mean)  
 library(shape)  
 Arrowhead(w_mean$lon, w_mean$lat, angle=alpha,   
      arr.length = 0.2, arr.type="curved")  
 lines(getMap(resolution = "low"), lwd=4)  


Wind connectivity between Azores and mainland


To complete this tutorial, we will use the downloaded wind data to compute wind connectivity between several points in mainland and Azores Archipelago. This kind of measures can be useful in order to establish theories about island colonization by plants or other organisms that could be influenced by wind currents. We will also use this connectivity to compute least cost paths between two locations, in order to virtually plot a wind highway connecting mainland with Azores islands.

 # First, we define our locations. 7 locations in the East   
 # and 7 in the West. The last one is a point in the Azores islands  
   
 loc <- matrix(c(-7, -4, -1, -9, -6,-10,-15,   
         -60, -56,-63,-76,-76,-81,-80, -27,  
         55, 50, 45, 40, 35, 30, 25,   
         55, 50, 45, 40, 35, 30, 25, 39), 15, 2)  
   
 colnames(loc) <- c("lon", "lat")  
 rownames(loc) <- c(seq(1,14,1), "Azores")  
   
 # Check how "loc" is looking   
   
 tail(loc)  
   
 # You can plot them  
 #plot(r_mean$wind.speed, col=acol(1000))  
 #lines(getMap(resolution = "low"), lwd=2)  
 #points(loc, col="red4", pch = 19, cex = 2)  


Now, we'll execute the next tasks:
  1. We'll use wind2raster to obtain raster layers of wind direction and speed of each day in the time series.
  2. Then, for each day, we'll compute the conductance matrix, a matrix with connectivity values between every single cell of the entire study area.
  3. Later, with each conductance matrix, we'll compute the cost between our locations and will store them in a cost_list object.
  4. Finally, for each day, if the Cost is not Infinite, we'll get the shortest path between two selected locations, and we'll store it into a “paths” object.

 # 1) Raster layers from wind time series  
   
 layers <- wind2raster(ww)  
   
 # 2) Conductance from wind direction and speed  
   
 Conductance <- flow.dispersion(layers, type="passive",  
                 output="transitionLayer")  
   
 # 3) and 4) Using a loop to compute cost between  
 # locations and least cost paths between location 11  
 # and Azores islands (as an example).  
   
 cost_list <- array(NA_real_, dim=c(15,15,30))  
 paths <- list(1:30)  
   
 library(gdistance)  
   
 for (i in 1:30){  
  cost_list[,,i] <- costDistance(Conductance[[i]],loc)  
  if (costDistance(Conductance[[i]],loc[11,] , loc[15,]) != Inf){  
   paths[[i]] <- shortestPath(Conductance[[i]], loc[11,], loc[15,],  
                 output="SpatialLines")  
   print(i)  
  }  
 }


Now, we will manage a little bit the results to plot them in a cool way.



 # Transform costs in connectivity  
   
 connectivity <- 1/cost_list  
   
 # Obtain connectivity averages  
   
 conn_avg <- apply(connectivity, c(1, 2), mean, na.rm = TRUE)  
   
 # Here we filter only connectivity from mainland to Azores  
 # As connectivity has arbitrary magnitude and units, we will   
 # multiply connectivity values 1000 times to avoid decimal zeros.  
   
 conn_azores <- conn_avg[,15]  
 conn_azores <- conn_azores *1000  
   
 # Finally, we can plot connectivity values relative to the width  
 # of the connections between mainland and Azores. You can power  
 # the values to increase the differences between them.  
   
 plot(r_mean$wind.speed, main= "Connectivity from mainland to Azores",
     col=acol(1000))  
 lines(getMap(resolution = "low"), lwd=1)  
 for (i in 1:14){  
  lines(rbind(loc[i,], loc[15,]), lwd = conn_azores[i]^3, col="red4")  
 }  
 points(loc, col="red", pch = 19, cex = 1.5)  
 text( loc+2, labels=c(round(conn_azores[1:14], 2),"Azores"),  
    col="white", cex= 1.5, font=2)  

Last but not least, we can use the shortest paths stored for each day to plot a density heatmap. This can help us to visualize wind corridors across our study area (our wind highways!).

 # We'll use spatstat library to convert the shortest  
 # paths lines into an psp object.  
   
 library(spatstat)  
   
 paths_clean <- paths[!sapply(paths, is.null)]  
 paths_merged <- paths_clean[[1]]  
 for (h in 2:length(paths_clean)) {  
  paths_merged <- rbind(paths_merged,paths_clean[[h]])  
 }  
 paths_psp <- as(paths_merged, "psp")  
   
 # Now, we can obtain a kernel density of the paths  
 # and convert it into a raster layer that can be plotted.  
 # You can play with sigma value to modify the smoothing.  
   
 lines_kernel <- density(paths_psp, sigma=1, dimyx=c(350,410))  
 kernel <- raster(lines_kernel)  
 #kernel[kernel<(maxValue(kernel)*0.1)]<-NA  
   
 ext <- extent(r_mean$wind.speed)  
 kernel <- extend(kernel, ext)  
 kernel[is.na(kernel)] <- minValue(kernel)  
   
 acol <- colorRampPalette(c("grey40","darkblue","red2",
                               "orange","yellow","white"))  
 plot(kernel, col=acol(1000), zlim=c(-0.01,5),  
       main = "Least cost paths density", xlab="Longitude",  
       ylab="Lattitude")  
 lines(getMap(resolution = "high"), lwd=2)





Isn't it cool?

 That's all for now!





7 comentarios:

  1. this is great. Thanks. (there is a mistake in the abstract of your rather cool looking paper "langauge" -> "language". Best, Mark

    ResponderEliminar
  2. Congrats, I really enjoy your exercise and package. Few details: there are some functions from raster package used, I think it will be useful to load at the beginning; colorRampPalette is also from grDevices package as it is used just once could be useful to load with ::
    Finally, the ylabel must be latitude instead of lattitude.
    Nice!!
    Walter

    ResponderEliminar
    Respuestas
    1. Thanks Walter!

      You are right, sorry about that! It's hard sometimes to fix all the code to be completely reproducible for all users... Feedback is essential!

      Cheers!

      Eliminar
    2. Javier, un placer haber contribuido con el trabajo, a las órdenes,
      Walter

      Eliminar
  3. Este comentario ha sido eliminado por el autor.

    ResponderEliminar
  4. image.plot(wind_layer$wind.speed, main="least cost paths by wind direction and speed",
    + col=terrain.colors(10), xlab="Longitude", ylab="Lattitude", zlim=c(0,7))
    Error in .local(x, ...) : invalid layer names
    ## Bonsoir depuis quelque temps je n'arrive plus à ploter ça pourquoi? chaque fois j'ai : Error in .local(x, ...) : invalid layer names comment resoudre ? Merci d'avance

    ResponderEliminar