R中的空间数据分析

> library(sp)
> library(maptools)
> library(raster)
> library(rgeos)

> maxd3 = readAsciiGrid("max.d3.txt")    # 加载数据
> raster .maxd3 = raster(maxd3)      # 转化为Raster
> maxd1 = readAsciiGrid("max.d1.txt")
> raster .maxd1 = raster(maxd1)     
> raster .max.d3.d1 = raster .maxd3 - raster .maxd1   # 比较3天与1天的数据
> rc <- reclassify(raster .max.d3.d1, c(-Inf,0,1, 0,Inf,NA))  # 对比较结果进行分类，提取出 有问题的数据

> raster .maxRange = raster .maxd1>0   #提取范围
> maxRange = rasterToPolygons(raster .maxRange, dissolve=TRUE)

> cols = grey.colors(20, 0.1, 0.9, 2.2)
> image(rc, col=cols, useRaster=TRUE)    #显示有问题的数据
> plot(maxRange, add=TRUE)               #显示整个范围