R代码求助

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

library(maps)
library(sp)
library(fields)
library(foreign)

"%&%"<-function(x,y)paste(x,y,sep="")

setwd("G:/2018Spring/econometrics/data/replication/data")

data <- read.table("yield_piecewise.txt",header=TRUE)
d=data[is.na(data$dday0C_smooth2000)==F & data$longitude>(-100) & is.na(data$longitude[i])==F,]

tchg = d$tavg_smooth2000 - d$tavg_smooth1980
tchg_sd = tchg/d$tavg_sd
pchg = (d$prec_smooth2000 - d$prec_smooth1980)/d$prec_smooth1980


pdf(file="G:/2018Spring/econometrics/data/replication/output/Figure1.pdf",height=6,width=12)
#temperature first
layout(matrix(c(1,1,1,2,3,3,3,4,5,5,5,6),nrow=4,byrow=F))
par(mar=c(0,0,1,1))
map(database="state",xlim=c(-100,-67))
title(main="Temperature",cex.main=2)
vals=round(tchg*100)
colz=tim.colors(max(vals,na.rm=T)-min(vals,na.rm=T))
vals = vals+abs(min(vals,na.rm=T)) #recenter so they are positive
for (i in 1:dim(data)[1]) {
#  if (is.na(data$longitude[i])==F & data$longitude[i]>(-100) ) {
    points(data$longitude[i],data$latitude[i],col=colz[vals[i]],pch=19,cex=1.3) #}
}
map(database="state",add=T)
vals=round(tchg,3)
mx=max(tchg,na.rm=T)my
mn=min(tchg,na.rm=T)
lng=(mx-mn)/302.54*
rg=seq(mn,mx,lng)
colz=tim.colors(length(rg))
par(mar=c(4,5,0,5))
hist(tchg,breaks=rg,col=colz,xlab="tobservations emp change (C)",main="",yaxt="n",ylab="",cex.axis=1.5,cex.lab=1.5)

#precip
toplot = pchg*100
toplot[toplot<(-50)]=-50
toplot[toplot>45]=45
map(database="state",xlim=c(-100,-67))
title(main="Precipitation",cex.main=2)
vals=round(toplot)
colz=rev(tim.colors(max(vals,na.rm=T)-min(vals,na.rm=T)))
vals = vals+abs(min(vals,na.rm=T)) #recenter so they are positive
for (i in 1:dim(data)[1]) {
  if (is.na(data$longitude[i])==F & data$longitude[i]>(-100) ) {
    points(data$longitude[i],data$latitude[i],col=colz[vals[i]],pch=19,cex=1.3)}
}
map(database="state",add=T)
vals=round(toplot,3)
mx=max(toplot,na.rm=T)
mn=min(toplot,na.rm=T)
lng=(mx-mn)/30
rg=seq(mn,mx,lng)
colz=rev(tim.colors(length(rg)))
par(mar=c(4,5,0,5))
hist(toplot,breaks=rg,col=colz,xlab="precip change (%)",main="",yaxt="n",ylab="",cex.axis=1.5,cex.lab=1.5)

#cornyield
ychg = log(data$cornyield_smooth2000) - log(data$cornyield_smooth1980)
map(database="state",xlim=c(-100,-67))
title(main="Corn yield",cex.main=2)
ychg[ychg<(-0.5)]=-0.5  #censor bottom at -0.5 to make plot prettier
vals=round(ychg*100)
colz=tim.colors(max(vals,na.rm=T)-min(vals,na.rm=T))
vals = vals+abs(min(vals,na.rm=T)) #recenter so they are positive
vals=max(vals,na.rm=T) - vals
for (i in 1:dim(data)[1]) {
  if (is.na(data$longitude[i])==F & data$longitude[i]>(-100) ) {
    points(data$longitude[i],data$latitude[i],col=colz[vals[i]],pch=19,cex=1.3)}
}
map(database="state",add=T)
vals=round(ychg,3)
mx=max(ychg,na.rm=T)
mn=min(ychg,na.rm=T)
lng=(mx-mn)/30
rg=seq(mn,mx,lng)
colz=rev(tim.colors(length(rg)))
par(mar=c(4,5,0,5))
hist(ychg,breaks=rg,col=colz,xlab="change in yield (log)",main="",yaxt="n",ylab="",cex.axis=1.5,cex.lab=1.5)

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：R代码 observations observation longitude database