> library(sp)
> library(maptools)
> library(raster)
> library(rgeos)
> maxd3 = readAsciiGrid("max.d3.txt") # 加载数据
> rasterlayer.maxd3 = raster(maxd3) # 转化为Raster
> maxd1 = readAsciiGrid("max.d1.txt")
> rasterlayer.maxd1 = raster(maxd1)
> rasterlayer.max.d3.d1 = rasterlayer.maxd3 - rasterlayer.maxd1 # 比较3天与1天的数据
> rc <- reclassify(rasterlayer.max.d3.d1, c(-Inf,0,1, 0,Inf,NA)) # 对比较结果进行分类，提取出 有问题的数据
> rasterlayer.maxRange = rasterlayer.maxd1>0 #提取范围
> maxRange = rasterToPolygons(rasterlayer.maxRange, dissolve=TRUE)
> cols = grey.colors(20, 0.1, 0.9, 2.2)
> image(rc, col=cols, useRaster=TRUE) #显示有问题的数据
> plot(maxRange, add=TRUE) #显示整个范围