如何绘制日食界限图？.doc

如何绘制日食界限图？ 张金龙 Jinlongzhang01@ 日食界限图是表示日食可见区域的一种地图。包括日食的南北界限，日出日落时刻线，日出日落时初亏线，日出日落时复圆线。如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) t