Revision fc70e7d5
Added by Adam Wilson over 10 years ago
| climate/research/cloud/MOD09/NDP-026D.R | ||
|---|---|---|
| 8 | 8 |
library(rasterVis) |
| 9 | 9 |
library(rgdal) |
| 10 | 10 |
library(reshape) |
| 11 |
|
|
| 11 |
library(maptools) |
|
| 12 |
library(rgeos) |
|
| 12 | 13 |
|
| 13 | 14 |
## Data available here http://cdiac.ornl.gov/epubs/ndp/ndp026d/ndp026d.html |
| 14 | 15 |
|
| ... | ... | |
| 51 | 52 |
cld=cld[,!grepl("Fq|AWP|NC",colnames(cld))]
|
| 52 | 53 |
cld$Amt[cld$Amt<0]=NA |
| 53 | 54 |
cld$Amt=cld$Amt/100 |
| 55 |
cld=cld[!is.na(cld$Amt),] |
|
| 56 |
|
|
| 57 |
## table of stations with > 20 observations per month |
|
| 58 |
cast(cld,StaID~YR,value="Nobs") |
|
| 59 |
mtab=ddply(cld,c('StaID','month'),function(df){ data.frame(count=sum(df$Nobs>20,na.rm=T))})
|
|
| 60 |
mtab2=mtab[ |
|
| 61 |
table(mtab$count>10) |
|
| 62 |
stem(mtab$count) |
|
| 54 | 63 |
|
| 55 | 64 |
## calculate means and sds for full record (1970-2009) |
| 56 | 65 |
Nobsthresh=20 #minimum number of observations to include |
| ... | ... | |
| 61 | 70 |
StaID=x$StaID[1], |
| 62 | 71 |
cld_all=mean(x$Amt[x$Nobs>=Nobsthresh],na.rm=T), # full record |
| 63 | 72 |
cldsd_all=sd(x$Amt[x$Nobs>=Nobsthresh],na.rm=T), |
| 73 |
cldn_all=length(x$Amt[x$Nobs>=Nobsthresh]), |
|
| 64 | 74 |
cld=mean(x$Amt[x$YR>=2000&x$Nobs>=Nobsthresh],na.rm=T), #only MODIS epoch |
| 65 |
cldsd=sd(x$Amt[x$YR>=2000&x$Nobs>=Nobsthresh],na.rm=T))})) |
|
| 66 |
cldm[,c("lat","lon")]=coordinates(st)[match(cldm$StaID,st$id),c("lat","lon")]
|
|
| 75 |
cldsd=sd(x$Amt[x$YR>=2000&x$Nobs>=Nobsthresh],na.rm=T), |
|
| 76 |
cldn=length(x$Amt[x$YR>=2000&x$Nobs>=Nobsthresh]))})) |
|
| 77 |
|
|
| 78 |
cldm[,c("lat","lon")]=coordinates(st)[match(cldm$StaID,st$id),c("lat","lon")]
|
|
| 67 | 79 |
|
| 68 | 80 |
|
| 69 | 81 |
|
| ... | ... | |
| 130 | 142 |
if(!file.exists("../teow/biomes.shp")){
|
| 131 | 143 |
teow=readOGR("/mnt/data/jetzlab/Data/environ/global/teow/official/","wwf_terr_ecos")
|
| 132 | 144 |
teow=teow[teow$BIOME<90,] |
| 133 |
biome=unionSpatialPolygons(teow,teow$BIOME, threshold=5)
|
|
| 145 |
biome=unionSpatialPolygons(teow,paste(teow$REALM,teow$BIOME,sep="_"), threshold=5)
|
|
| 134 | 146 |
biomeid=read.csv("/mnt/data/jetzlab/Data/environ/global/teow/official/biome.csv",stringsAsFactors=F)
|
| 135 |
biome=SpatialPolygonsDataFrame(biome,data=biomeid[as.numeric(row.names(biome)),]) |
|
| 147 |
realmid=read.csv("/mnt/data/jetzlab/Data/environ/global/teow/official/realm.csv",stringsAsFactors=F,na.strings = "TTTT")
|
|
| 148 |
dt=data.frame(code=row.names(biome),stringsAsFactors=F) |
|
| 149 |
dt[,c("realmid","biomeid")]=do.call(rbind,strsplit(sub(" ","",dt$code),"_"))
|
|
| 150 |
dt$realm=realmid$realm[match(dt$realmid,realmid$realmid)] |
|
| 151 |
dt$biome=biomeid$BiomeID[match(dt$biomeid,biomeid$Biome)] |
|
| 152 |
row.names(dt)=row.names(biome) |
|
| 153 |
biome=SpatialPolygonsDataFrame(biome,data=dt) |
|
| 136 | 154 |
writeOGR(biome,"../teow","biomes",driver="ESRI Shapefile",overwrite=T) |
| 137 | 155 |
} |
| 138 | 156 |
biome=readOGR("../teow/","biomes")
|
| 139 | 157 |
projection(biome)=projection(st) |
| 140 | 158 |
#st$biome=over(st,biome,returnList=F)$BiomeID |
| 141 | 159 |
dists=apply(gDistance(st,biome,byid=T),2,which.min) |
| 142 |
st$biome=biome$BiomeID[dists] |
|
| 160 |
st$biomec=biome$code[dists] |
|
| 161 |
st$realm=biome$realm[dists] |
|
| 162 |
st$biome=biome$biome[dists] |
|
| 143 | 163 |
|
| 164 |
cldm$biomec=st$biomec[match(cldm$StaID,st$id)] |
|
| 165 |
cldm$realm=st$relam[match(cldm$StaID,st$id)] |
|
| 144 | 166 |
cldm$biome=st$biome[match(cldm$StaID,st$id)] |
| 145 | 167 |
|
| 146 | 168 |
|
Also available in: Unified diff
updating figures