############################################################# #####GET SISTER PAIR RANGE SIZE AND CO_OCCURENCE SUMMARY ### ############################################################# ###################### # step 1. First run sisters script to get list of sisters (sister.stats.function) and then import summary table ###################### sis=read.csv(file ="summary.table.csv" , as.is=TRUE, sep=",") rownames(sis)=sis$sister ###################### ### step 2. save function at bottom of script to memory. ###################### ###################### ##### step 3. get list of xy, yx sister taxa######### ###################### tmp=strsplit(rownames(sis),split=" ") tmp2 = do.call(rbind,tmp) sis$ab = paste(tmp2[,1],tmp2[,2],sep="_") sis$ba = paste(tmp2[,2],tmp2[,1],sep="_") head(sis) library(raster) ###################### ##### step 4. get geography data for sister taxa######### ###################### ### note: this is super slow amsinckia=sis.geo("amsinckia") arabidopsis =sis.geo("arabidopsis") capsella =sis.geo("capsella") clarkia =sis.geo("clarkia") collinsia =sis.geo("collinsia") downingia =sis.geo("downingia") erodium=sis.geo("erodium") #no gbif records for gaillardotii, lucidum, cossonii, pelargoniflorum, dalechampia =sis.geo("dalechampia") lasthenia =sis.geo("lasthenia") leptosiphon =sis.geo("leptosiphon") limnanthes =sis.geo("limnanthes") medicago =sis.geo("medicago") mimulus =sis.geo("mimulus") oenothera.anogra =sis.geo("oenothera.anogra") oenothera.oenothera =sis.geo("oenothera.oenothera") polemonium =sis.geo("polemonium") primula =sis.geo("primula") # #no gbif records for alcalina, frondosa saltugilia =sis.geo("saltugilia") schiedea =sis.geo("schiedea") # no gbif records for attenuata, kaalae (lat was wrong), mannii, verticillata schizanthus =sis.geo("schizanthus") # no gbif records for porrigens, parvulus, tricolor dat=rbind(amsinckia,arabidopsis,capsella,clarkia,collinsia,dalechampia,downingia,erodium,lasthenia,leptosiphon,limnanthes,mimulus,medicago,oenothera.anogra,oenothera.oenothera,polemonium,primula, saltugilia, schiedea, schizanthus) colnames(dat) tmp = subset(dat, select = -c(sp.A.y,sp.B.y)) colnames(tmp) colnames(tmp)[colnames(tmp)=="sp.A.x"] <- "sp.A" colnames(tmp)[colnames(tmp)=="sp.B.x"] <- "sp.B" #write.table(tmp,"summary.table.csv",sep=",",row.names=FALSE,col.names=TRUE) ####################################### #####function to get geography data for sister taxa ######### ####################################### sis.geo = function(GENUS){ # Make a vector of all pairwise comparisons between species IndPWComps = function(our.inds){ ind.comps = do.call(rbind,lapply(1:(length(our.inds)-1),function(X){ cbind(our.inds[X],our.inds[(X+1):length(our.inds)]) })) comp.list =lapply(1:nrow(ind.comps),function(X){ind.comps[X,]}) names(comp.list) = paste(ind.comps[,1],ind.comps[,2],sep="_") return(comp.list) } # Import a dataframe for our genus GenusData = function(genus.filename){ temp.data = read.csv(file = genus.filename, as.is=TRUE, sep=",") return(temp.data) } # Get the area of cells occupied by each species and jointly occupied by both species for all species combinations GetRangesSummary = function(genus.data,res.rx){ rx <- raster() #create a raster layer res(rx) <- res.rx #set it's resolution at 1 km sp.area = do.call(c, lapply(unique(genus.data$Species),function(sp){ my.xy=genus.data[genus.data$Species == sp,c("Longitude","Latitude")] #get lat longs for species my.r <- rasterize(my.xy, rx, field=1) #set raster cells where species present = 1 my.r.area=raster::area(my.r,na.rm=TRUE) #makes a raster with area values for cells occupied by sp sum(as.matrix(my.r.area),na.rm=TRUE) #get area sum of cells occupied by sp })) names(sp.area) = unique(genus.data$Species) sp.pairwise.list = IndPWComps(unique(genus.data$Species)) genus.occupancy = lapply(unique(genus.data$Species),function(sp){ unique(cellFromXY(rx,genus.data[genus.data$Species == sp,c("Longitude","Latitude")]))}) names(genus.occupancy) = unique(genus.data$Species) sp.overlap = do.call(rbind, lapply(sp.pairwise.list,function(pw.sps){ A.area = sp.area[[pw.sps[1]]] B.area= sp.area[[pw.sps[2]]] A.geo = genus.occupancy[[pw.sps[1]]] #cells with spA B.geo = genus.occupancy[[pw.sps[2]]] #cells with spB overlap = intersect (A.geo, B.geo) #cells with both species net.geo = unique(append(A.geo,B.geo)) #cells with spA, spB, or both A.coords = subset(genus.occur, Species==pw.sps[1]) B.coords = subset(genus.occur, Species==pw.sps[2]) my.xy=xyFromCell(rx, overlap) #get lat longs for overlap if (nrow(my.xy)<1) { total.overlap=0 } else { my.r <- rasterize(my.xy, rx, field=1) #set raster cells where species present = 1 my.r.area=raster::area(my.r,na.rm=TRUE) #makes a raster with area values for cells occupied by both species total.overlap=sum(as.matrix(my.r.area),na.rm=TRUE) #get area sum of cells occupied by both species } net.xy=xyFromCell(rx, net.geo) #get lat longs for net area if (nrow(net.xy)<1) { net.area=0 } else { net.r <- rasterize(net.xy, rx, field=1) #set raster cells where species present = 1 net.r.area=raster::area(net.r,na.rm=TRUE) #makes a raster with area values for cells occupied by both species net.area=sum(as.matrix(net.r.area),na.rm=TRUE) #get area sum of cells occupied by both species } require(fields) distMatrix <- rdist.earth(A.coords[,c('Longitude','Latitude')], B.coords[,c('Longitude','Latitude')], miles = TRUE) min.dist = min(distMatrix) summary=c(A.range = A.area, B.range = B.area, total.overlap = total.overlap, net.area=net.area, min.dist = min.dist) }) ) return(sp.overlap) } # Get range and range overlap across resolutions GetRangeSummaries=function(genus.data,res.rxs){ these.summaries = lapply(res.rxs,function(res.rx){GetRangesSummary(genus.data,res.rx)}) names(these.summaries) = paste("res",round(res.rxs,digits=3),sep=":") return(these.summaries) } #After entering functions to memory run these scripts genus.occur = GenusData(sprintf("data geographic/ft.clean.gbif.%s.csv",GENUS)) #our file, with dates from 0-365 genus.area.summary = GetRangeSummaries(genus.occur,c(1,0.5,0.1,0.05)) #A list with three things: #make overlap stats min.dist=genus.area.summary$'res:1'[,5] A.range.res1=genus.area.summary$'res:1'[,1] B.range.res1=genus.area.summary$'res:1'[,2] overlap.net.res1=genus.area.summary$'res:1'[,3]/(genus.area.summary$'res:1'[,4]) overlap.abs.res1=genus.area.summary$'res:1'[,3] overlap.A.res1=genus.area.summary$'res:1'[,3]/genus.area.summary$'res:1'[,1] overlap.B.res1=genus.area.summary$'res:1'[,3]/genus.area.summary$'res:1'[,2] A.range.res0.5=genus.area.summary$'res:0.5'[,1] B.range.res0.5=genus.area.summary$'res:0.5'[,2] overlap.net.res0.5=genus.area.summary$'res:0.5'[,3]/(genus.area.summary$'res:0.5'[,4]) overlap.abs.res0.5=genus.area.summary$'res:0.5'[,3] overlap.A.res0.5=genus.area.summary$'res:0.5'[,3]/genus.area.summary$'res:0.5'[,1] overlap.B.res0.5=genus.area.summary$'res:0.5'[,3]/genus.area.summary$'res:0.5'[,2] A.range.res0.1=genus.area.summary$'res:0.1'[,1] B.range.res0.1=genus.area.summary$'res:0.1'[,2] overlap.net.res0.1=genus.area.summary$'res:0.1'[,3]/(genus.area.summary$'res:0.1'[,4]) overlap.abs.res0.1=genus.area.summary$'res:0.1'[,3] overlap.A.res0.1=genus.area.summary$'res:0.1'[,3]/genus.area.summary$'res:0.1'[,1] overlap.B.res0.1=genus.area.summary$'res:0.1'[,3]/genus.area.summary$'res:0.1'[,2] A.range.res0.05=genus.area.summary$'res:0.05'[,1] B.range.res0.05=genus.area.summary$'res:0.05'[,2] overlap.net.res0.05=genus.area.summary$'res:0.05'[,3]/(genus.area.summary$'res:0.05'[,4]) overlap.abs.res0.05=genus.area.summary$'res:0.05'[,3] overlap.A.res0.05=genus.area.summary$'res:0.05'[,3]/genus.area.summary$'res:0.05'[,1] overlap.B.res0.05=genus.area.summary$'res:0.05'[,3]/genus.area.summary$'res:0.05'[,2] #split pairwise names into two columnes sp.PW <-data.frame(do.call('rbind',strsplit(row.names(genus.area.summary$'res:1'),'_',fixed=TRUE))) sp.PW=subset(sp.PW,select=c(X1,X2)) colnames(sp.PW)=c("sp.A","sp.B") #put overlap stats and mating diffs into table dat=as.data.frame(cbind(sp.PW,min.dist,A.range.res1,B.range.res1,A.range.res0.5,B.range.res0.5,A.range.res0.1,B.range.res0.1,A.range.res0.05,B.range.res0.05,overlap.net.res1,overlap.abs.res1,overlap.A.res1,overlap.B.res1,overlap.net.res0.5,overlap.abs.res0.5,overlap.A.res0.5,overlap.B.res0.5,overlap.net.res0.1,overlap.abs.res0.1,overlap.A.res0.1,overlap.B.res0.1,overlap.net.res0.05,overlap.abs.res0.05,overlap.A.res0.05,overlap.B.res0.05)) #select sister pairs tmp=subset(sis,genus==sprintf("%s",GENUS)) #make xy and yx columns of pair names on geography dat$ab=paste(dat[,1], dat[,2],sep="_") #select and merge geography of just sister pairs t1=merge(tmp, dat,by="ab") t2=merge(tmp, dat,by.x="ba",by.y="ab") sis.geo=rbind(t1,t2) }