library(ggmap) library(mapdata) # eqdata<-read.table(text="",col.names=c("time", "longitude", "latitude", "depth", "mag")) ymd<-seq(as.Date("2016-04-14"),as.Date("2016-05-02"), by="1 day") for (i in 1:length(ymd)){ date<-gsub("-","",ymd[i]) #ウェブスクレイピング url<-paste0("http://www.data.jma.go.jp/svd/eqev/data/daily_map/",date,".html") # doc<-readLines(url,encoding ="UTF-8") kensaku<-paste0(substr(date,1,4)," ",gsub("\\<0"," ",substr(date,5,6))," ",gsub("\\<0"," ",substr(date,7,8))) # x<-doc[grep(paste0("\\<",kensaku),doc)] #全角文字 ° が入っているので半角スペース２つに変換する x<-gsub("°", " ",x) time<-paste0(substr(x,1,4),"-",substr(x,6,7),"-",substr(x,9,10)," ",substr(x,12,16),":",substr(x,18,21)) # #南緯、西経の場合も考慮に入れると latitude = as.numeric(paste0(substr(x,23,23),as.character(as.numeric(substr(x,24,25))+as.numeric(substr(x,28,29))/60+as.numeric(substr(x,31,31))/600))) longitude = as.numeric(paste0(substr(x,34,34),as.character(as.numeric(substr(x,35,37))+as.numeric(substr(x,40,41))/60+as.numeric(substr(x,43,43))/600))) depth<-as.numeric(substr(x,48,50)) mag<-as.numeric(substr(x,55,58)) eq<-data.frame(time,longitude,latitude,depth,mag) eqdata<-rbind(eqdata,subset(eq,mag>=-4)) } # #重複した行を取り除く eqdata<-unique(eqdata) #head(eqdata) ; tail(eqdata) #save(eqdata,file="eqdata20160414_0502.RData") #load("eqdata20160414_0502.RData") # #震源の深さ # depth<10km eqdata1<-subset(eqdata,depth<10) eqdata1$depth_rank<-" 10km 未満" eqdata1$col<-"red" # 10<=depth<20 eqdata2<-subset(eqdata,depth>=10 & depth<20) eqdata2$depth_rank<-" 10km以上20km未満" eqdata2$col<-"orange" # 20<=depth<40 eqdata3<-subset(eqdata,depth>=20 & depth<40) eqdata3$depth_rank<-" 20km以上40km未満" eqdata3$col<-"gold4" # 40<=depth<80 eqdata4<-subset(eqdata,depth>=40 & depth<80) eqdata4$depth_rank<-" 40km以上80km未満" eqdata4$col<-"green" # 80<=depth<150 eqdata5<-subset(eqdata,depth>=80 & depth<150) eqdata5$depth_rank<-" 80km以上150km未満" eqdata5$col<-"blue" # 150<=depth eqdata6<-subset(eqdata,depth>=150) eqdata6$depth_rank<-"150km以上" eqdata6$col<-"purple" # eq<-rbind(eqdata1,eqdata2,eqdata3,eqdata4,eqdata5,eqdata6) #並べ替え:マグニチュードの昇順 sortlist <- order(eq$mag) eq<- eq[sortlist,] ## #mag>=1の地震を抽出 subeq<-subset(eq,mag>=1) #mag<1 mag1und<-subset(eq,mag<1) #最大マグニチュードの地震の震源を地図の中心にする。 subset(eqdata,mag==max(eqdata$mag)) # # time longitude latitude depth mag #5312 2016- 4-16 01:25: 5.4 130.7617 32.75333 12 7.3 # par(family = "TakaoExMincho") #png("eq20160414_0502_01.png",width=1000,height=800) #最大マグニチュードの地震の震源を地図の中心にする場合（１つしか該当する地震がないとき） gmap <- get_map(c(lon = subset(eqdata,mag==max(eqdata$mag))$longitude,subset(eqdata,mag==max(eqdata$mag))$latitude),zoom = 9, maptype = "terrain") #地図の中心を指定する場合 #gmap <- get_map(c(130.7617,32.75333),zoom = 9, maptype = "terrain") attr(gmap,"b") # ll.lat ll.lon ur.lat ur.lon #1 32.00995 129.8841 33.48828 131.6419 # p<-ggmap(gmap)+ #mag>=1の地震をプロット geom_point(data=subeq, aes(x=longitude, y=latitude,fill=depth_rank,size=mag),alpha=0.8,shape=21,color="gray20")+ #プロットする色を指定したい場合 scale_fill_manual(values=c("red","orange","gold4","green","blue","purple"), breaks = c(" 10km 未満"," 10km以上20km未満"," 20km以上40km未満"," 40km以上80km未満"," 80km以上150km未満","150km以上"))+ #プロットする大きさの範囲を指定したい場合 scale_size(breaks=c(1,3,5,7),range = c(1,10)) + labs(title=paste0(ymd[1],"～",ymd[length(ymd)]," 地震 (mag>=1)")) # #png("eq20160414_0502_01.png",width=1000,height=800) p #dev.off()
| 
