library(raster) source("http://www.fttsus.jp/worldmesh/R/worldmesh.R") # Gen_frame_e_to_w<-function(efgs_code){ Y <- as.numeric(gsub('.*N([0-9]+)[EW].*', '\\1', efgs_code))*1000 X <- as.numeric(gsub('.*[EW]([0-9]+)', '\\1', efgs_code))*1000 frame <- as(extent(X,X+1000,Y,Y+1000), "SpatialPolygons") proj4string(frame) <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs") frame_etow <- spTransform(frame,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")) return(frame_etow) } # Cal_cover_efgs_to_world<-function(efgs_code){ frame_etow<-Gen_frame_e_to_w(efgs_code) min_long<-xmin(frame_etow) min_lat<-ymin(frame_etow) max_long<-xmax(frame_etow) max_lat<-ymax(frame_etow) mid_lat<-(min_lat+max_lat)/2 mid_long<-(min_long+max_long)/2 NE=c(max_lat,max_long) NW=c(max_lat,min_long) SE=c(min_lat,max_long) SW=c(min_lat,min_long) M1=c(min_lat,mid_long) M2=c(max_lat,mid_long) M3=c(mid_lat,min_long) M4=c(mid_lat,max_long) check_points<-rbind(NE,NW,SE,SW,M1,M2,M3,M4) cover<-apply(check_points,1,function(p){return(cal_meshcode3(p[1],p[2]))}) return(unique(cover)) } # Plot_cover_and_efgs_code<-function(efgs_code){ frame_etow <-Gen_frame_e_to_w(efgs_code) cover<-Cal_cover_efgs_to_world(efgs_code) cover_grids<-lapply(cover,meshcode_to_latlong_grid) cover_frames<-lapply(cover_grids, function(cover){ frame<-as(extent(cover$long0,cover$long1,cover$lat1,cover$lat0), "SpatialPolygons") proj4string(frame) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") return(frame) } ) xmin<-xmin(frame_etow) xmax<-xmax(frame_etow) ymin<-ymin(frame_etow) ymax<-ymax(frame_etow) xd<-xmax-xmin yd<-ymax-ymin plot(frame_etow,col="blue",xlim=c(xmin-xd,xmax+xd),ylim=c(ymin-yd,ymax+yd)) print(cover) print(cover_frames) print(frame_etow) t<-lapply(cover_frames,function(x){plot(x,add=TRUE)}) } # efgs_code <- "1kmN3211E4322" Plot_cover_and_efgs_code(efgs_code)