Mapping

The goal here is to find overlap of highland (>1700m) and area under maize cultivation. Data are from worldclim and earthstat Load pacakges

require(sp)

## Loading required package: sp

## Warning: package 'sp' was built under R version 3.1.2

require(rgdal)

## Loading required package: rgdal
## rgdal: version: 0.9-1, (SVN revision 518)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.1, released 2014/09/24
## Path to GDAL shared files: /usr/local/share/gdal
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: (autodetected)

require(raster)

## Loading required package: raster

## Warning: package 'raster' was built under R version 3.1.1

require(ncdf)

## Loading required package: ncdf

## Warning: package 'ncdf' was built under R version 3.1.1

This plots the area in the highlands (yellow) and pixels (each pixel is a 10km x 10km square) in which >1% of the land is cultivated currently for maize.

#load data, plot area of interest in green
Altitude=raster("~/Desktop/area/alt.bil") #altitude data from worldclim
Maize=raster("~/Desktop/area/maize_AreaYieldProduction.nc",level=1) #maize cultivation data from earthstat; first level is percent land cultivated
e<-extent(c(-89.82036,-28.74251,-56.5066,12.90124)) #defines extent of South America
SAAlt<-crop(Altitude,e) #apply extent to altitude layer
SAMaize<-crop(Maize,e) #apply extent to maize layer
a<-SAAlt>1700 #define threshold for altitude
m<-SAMaize>0.01 #define threshold for proportion of pixel under maize cultivation
intersect<-a+m #add layers with raster algebra; layers meeting both criteria (altitude and maize cultivation have value of "2")
plot(intersect, main="Intersect of Maize Cultivation > 0 and Altitude >1700m")

Now we want to know in that area, exactly how many ha of maize fields are there.

#calculation of total area of grids with maize above 1700m
ivals<-values(intersect) #creates an object with just the values from the intersect
length(subset(ivals,ivals==2))*100 #gives the area in km^2 of intersect assuming 5 arc-minute resolution--100km^2/pixel

## [1] 270200

#calculation of extent of area currently under maize cultivation in grids
projection(SAAlt)<-"+proj=utm +zone=48 +datum=WGS84" #to stack layer each must be in the same geographic projection
projection(SAMaize)<-"+proj=utm +zone=48 +datum=WGS84" #to stack layer each must be in the same geographic projection
ST<-stack(SAAlt,SAMaize) #stacking altitude and maize layers so attributes can be retained in a query
STM<-as.matrix(ST) #converting S4 object to matrix
STD<-data.frame(STM) #converting matrix to data.frame for subset query
sub<-subset(STD, alt > 1700 & maizeData > 0.01) # identifying subset of grids with maize cultivation and >1700m
Akm2=sum(sub$maizeData)*100 #calculating area in km^2
A=Akm2*100 #ha

We now jump to Peter’s code. Waiting time to parallel mutation is just the inverse of probability of new mutation, which Peter shows is: \(\frac{2 \mu \rho A s_b}{\xi^2}\) where A is area, \(\rho\) planting density, and \(\xi^2\) the variance in reproductive success.

Here I am assuming a selection coefficient of \(10^-5\) for each mutation (weak selection on mutations underlying a quantitative trait), a mutation rate of \(mu=3 x 10^-8\) from Clark, and \(\xi^2=30\) (if we use 20 which is Peter’s estimate for wildtype none of the qualitative conclusions change).

rho=20000
sb=10^-5
xisq=30
mu=3*10^-8
Tmut=1/( mu * (2 * rho * A * sb)/xisq )
print(Tmut)

## [1] 4162.854