R中的空间数据分析

> 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)               #显示整个范围