2010-11-16

如何绘制日食界限图？

日食界限图是表示日食可见区域的一种地图。包括日食的南北界限，日出日落时刻线，日出日落时初亏线，日出日落时复圆线。如2011年1月4日将发生一次日偏食，在欧洲大部分，非洲撒哈拉及其以北地区，中亚，西亚以及我国的西部可以见到日偏食。可见区域由若干条曲线围成。各曲线代表的含义分别为，

I，日出时复圆。
II，日出时食甚
III，日出时初亏
IV，日落时复圆
V，日落时食甚
VI，日落时初亏
还有一条线段表示日食南界。在这个区域内部，可以见到日偏食的完整过程。区域之外则不可见。由于月亮自西向东运行，其速度远远超过太阳的向东运行的速度。因此，日食总是从西侧开始，在东侧结束。
公元前2000年到公元3000年每一次日食kmz文件，可以在 http://xjubier.free.fr/en/site_pages/SolarEclipsesGoogleEarth.html 下载。kmz文件可以被google earth读取。
为了生成清晰度更高的日食界限图，笔者编写了相应的R函数，以绘制日食界限图。该函数基于maptools、mapproj、mapdata。使用前请先下载这三个程序包。
由于绘制地图必须采用一定的投影方式，这里默认的投影方式为斜轴等面积投影。当然，高手可以通过阅读R的帮助信息，改变相应的地图投影方式，线段的颜色等，以达到满意的效果。
图2 - 6
plotkml绘制的日食界限图

plotkml <- 
function(file, projection = "azequalarea" , parameter = NULL, grid = TRUE, info = TRUE,... ){
	 require(maptools)
     require(mapproj)
     require(mapdata)
    #match.arg(projection)
    zz <- getKMLcoordinates(file, ignoreAltitude = TRUE)
	x11(11,9)
    m <- map('world',proj= projection, parameter = parameter, orient=c(zz[[1]][2], zz[[1]][1], 0), ...)
    
	for(i in 1:length(zz)){                              ## Add lines to graph
	    if(is.data.frame(zz[[i]])){
		    if(ncol(zz[[i]]) > 3){
		        i = i + 1
            }		                                     ## Avoid  ncol > 3
		}
		
		if(is.list(zz[[i]])){
		    i = i + 1			                         ## Avoid  list > 3
		}
		
        if(is.vector(zz[[i]])){
            points(mapproject(list(y = zz[[i]][2], x= zz[[i]][1])), pch = "*", cex = 2, col = 2)
		}
		
		if(is.matrix(zz[[i]])){
            lines(mapproject(list(y = zz[[i]][,2], x= zz[[i]][,1])), col = "blue", ...)
        } 
		## print(i) used in debugging
    }
    
	if(grid){   ## Add map grid
        map.grid(m, lty = 1, col = "grey")
    }
	file1 <- readLines(file)
	
	if(title){ ## Add title
	    #file <- "c:/ASE_1800_04_24__4cdcf4f5330fb.kml"
        daterecord <- file1[grepl("when", file1)]
        date <- substring(gsub(" ", "", daterecord), 7, 16)
        type <- file1[grepl("Type:", file1)]
        stype <- gsub(" ", "", type)
        typesolar <- substring(stype, 12, regexpr("<br", stype) - 1)
        title(paste(typesolar, "solar eclipse on", date))
	}
	if(info){  ## Add details for this eclipse
	    ##UT
		## detail information
		bg <- which(grepl("<li>Magnitude:", file1))
		info <- file1[bg:(bg+4)]
		infosub <- gsub("</li>", "", gsub("&nbsp;","",gsub("<li>", "",gsub(" ", "", info))))
        text(-1.3, seq(1.4, 1.1, by = -0.1),infosub[1:4])
		
		## Partial Eclipse Begins
		bg <- min(which(grepl("<name>P1", file1)))
        sub <- file1[bg:(bg+100)]
        ind <- min(which(regexpr("<longitude>",sub)>0))
		info <- sub[c(ind, ind +1)]
		info1 <- gsub("de", "de:",gsub("/longitude","",gsub("/latitude", "", gsub("[<->]", "", gsub(" ", "",info)))))
		information <- c("Partial Eclipse begins:", info1)
		text(-1.3,seq(-1.1, -1.3, by = -0.1),information)
		
		## Partial Eclipse Ends
		bg <- min(which(grepl("last external", file1)))
        sub <- file1[bg:(bg+100)]
        ind <- min(which(regexpr("<longitude>",sub)>0))
		info <- sub[c(ind, ind +1)]
		info1 <- gsub("de", "de:",gsub("/longitude","",gsub("/latitude", "", gsub("[<->]", "", gsub(" ", "",info)))))
		information <- c("Partial Eclipse ends:", info1)
		text(1.3,seq(-1.1, -1.3, by = -0.1),information)
		
		## Greatest Eclipse Point
		bg <- min(which(grepl("Greatest Eclipse Point", file1)))
        sub <- file1[bg:(bg+100)]
        ind <- min(which(regexpr("<longitude>",sub)>0))
		info <- sub[c(ind, ind +1)]
		info1 <- gsub("de", "de:",gsub("/longitude","",gsub("/latitude", "", gsub("[<->]", "", gsub(" ", "",info)))))
		information <- c("Greatest Eclipse:", info1)
		text(1.3,seq(1.3, 1.1, by = -0.1),information)
	}
}