R语言McSpatial_R语言：地图叠加与分析

setwd("d:/overlay")

library(maptools)

library(PBSmapping)

library(rgdal)

SouthAmBase = readShapePoly("SouthAmericaStates",proj4string=CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))

class(SouthAmBase)

plot(SouthAmBase)

print("here is the South America base map.....Hit key to continue")

SouthAmBaseProj = spTransform(SouthAmBase,CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))

SouthAmBaseProjPS = SpatialPolygons2PolySet(SouthAmBaseProj)

plotPolys(SouthAmBaseProjPS,proj = TRUE,col="wheat1",xlab="longitude",ylab="latitude")

BrazilProj   = SouthAmBaseProj[SouthAmBaseProj$COUNTRY == "BRASIL",]

BrazilProjPS = SpatialPolygons2PolySet(BrazilProj)

par(lwd = 2) # line width of 2, thicker than default of 1

addPolys(BrazilProjPS ,border="green")

attr(BrazilProjPS, "projection")

AreaBrazil = calcArea(BrazilProjPS,rollup=1)

BrazilArea = sum(AreaBrazil[1:length(AreaBrazil[,1]),2])

print(sprintf("Brazil area (green outline): %g sq km. Hit key to continue",BrazilArea))

plotPolys(SouthAmBaseProjPS,proj = TRUE,col="wheat1",xlab="longitude",ylab="latitude")

HabitatL1A = readShapePoly("SouthernArmadilloPoly",proj4string=CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))

HabitatL1Proj = spTransform(HabitatL1A,CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))

HabitatL1ProjPS = SpatialPolygons2PolySet(HabitatL1Proj)

print("Press  to view Southern Armadillo geographic range ..")

addPolys(HabitatL1ProjPS,col="red")

HabitatL2 = readShapePoly("GreaterArmadilloPoly",proj4string=CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))

HabitatL2Proj = spTransform(HabitatL2,CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))

HabitatL2ProjPS = SpatialPolygons2PolySet(HabitatL2Proj)

print("Press  to view Greater Armadillo geographic range  ..")

addPolys(HabitatL2ProjPS,col="green4")

HabitatL3 = readShapePoly("VesperMousePoly",proj4string=CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))

HabitatL3Proj = spTransform(HabitatL3,CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))

HabitatL3ProjPS= SpatialPolygons2PolySet(HabitatL3Proj)

print("Press  to view Vesper Mouse geographic range  ..")

addPolys(HabitatL3ProjPS,col="blue1")

HabitatL2B = readShapePoints("GreaterArmadilloPts",proj4string=CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))

HabitatL2PointProj = spTransform(HabitatL2B,CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))

nSightings = length(HabitatL2PointProj@data[,1])

SpeciesSightingEventSet

print(sprintf("There are %d recorded Greater Armadillo sightings.",nSightings ))

print("Press  to view Species sighting points..")

addPoints(SpeciesSightingEventSet,col="orange",pch=19)

用数据可以回答以下几个问题：

1) 每个种的分布区面积是多少?

#计算SouthernArmadillo分布区面积

attr(HabitatL1ProjPS, "projection")

AreaHab1 = calcArea(HabitatL1ProjPS,rollup=1)

SouthernArmadilloHabArea= sum(AreaHab1[1:length(AreaHab1[,1]),2])

print(sprintf("Southern Naked-Tailed Armadillo geographic range  area (Red): %g sq km",SouthernArmadilloHabArea))

#计算GreaterArmadillo分布区面积

attr(HabitatL2ProjPS, "projection")

AreaHab2 = calcArea(HabitatL2ProjPS,rollup=1)

GreaterArmadilloHabArea = sum(AreaHab2[1:length(AreaHab2[,1]),2])

print(sprintf("Greater Naked-Tailed Armadillo geographic range  area (blue): %g sq km",GreaterArmadilloHabArea))

#计算VesperMouse分布区面积

attr(HabitatL3ProjPS, "projection")

AreaHab3 = calcArea(HabitatL3ProjPS,rollup=1)

VesperMouseHabArea = sum(AreaHab3[1:length(AreaHab3[,1]),2])

print(sprintf("Vesper Mouse geographic range  area (Green): %g sq km",VesperMouseHabArea))

print(sprintf("There are %d recorded Greater Armadillo sightings (Orange).",nSightings))

print("Press  to demonstrate polygon overlay and analysis..")

问题二 三个种覆盖的总面积是多少？

print("Redrawing base map in new window..")

dev.set(1)

plotPolys(SouthAmBaseProjPS,proj = TRUE,col="wheat1",xlab="longitude",ylab="latitude")

print("Calculating the UNION of Armadillo Ranges (BLUE area)..")

Join1 = joinPolys(HabitatL1ProjPS,HabitatL2ProjPS,"UNION",maxVert = 1.0e+06)

print("...Here is the UNION of Armadillo Ranges..")

addPolys(Join1,col="blue")

print("Calculating the UNION of ALL THREE Ranges..")

AllHabitatUnion = joinPolys(Join1,HabitatL3ProjPS,"UNION",maxVert = 1.0e+06)

AllHabUnionCB = combinePolys(AllHabitatUnion)

print("Here is UNION of ALL THREE Ranges (RED area)..")

addPolys(AllHabUnionCB,col="red")

attr(AllHabUnionCB, "projection")

AreaHabitatUnion = calcArea(AllHabUnionCB,rollup=1)

AllHabitatArea = sum(AreaHabitatUnion[1:length(AreaHabitatUnion[,1]),2])

print(sprintf("Area covered by one or more geographic range: %g sq km",AllHabitatArea))

问题三 三个种同时存在的面积是多少？

dev.set(1)

plotPolys(SouthAmBaseProjPS,proj = TRUE,col="wheat1",xlab="longitude",ylab="latitude")

AllTwoHabitatArea   = 0.0

AllThreeHabitatArea = 0.0

print("Calculating the INTERSERCTION of FIRST Two Ranges..")

Join1 = joinPolys(HabitatL1ProjPS,HabitatL2ProjPS,"INT",maxVert = 1.0e+06)

if (length(Join1) > 0)

{

HabL1CB = combinePolys(Join1)

attr(HabL1CB, "projection")

AreaHabInt = calcArea(HabL1CB,rollup=1)

AllTwoHabitatArea = sum(AreaHabInt[1:length(AreaHabInt[,1]),2])

addPolys(HabL1CB ,col="red")

print("Calculating the INTERSERCTION of LAST Two Ranges..")

Join3 = joinPolys(HabL1CB,HabitatL3ProjPS,"INT",maxVert = 1.0e+06)

if (length(Join3) > 0)

{

HabL3CB = combinePolys(Join3)

attr(HabL3CB, "projection")

AreaHabInt = calcArea(HabL3CB,rollup=1)

AllFouroHabitatArea = sum(AreaHabInt[1:length(AreaHabInt[,1]),2])

addPolys(HabL3CB,col="blue")

}

}

print(sprintf("Area (RED)   covered by the first two geographic range: %g sq km",AllTwoHabitatArea))

print(sprintf("Area (GREEN) covered by the all three geographic range: %g sq km",AllThreeHabitatArea))

问题四

问题四 巴西有多大的面积被某个种的所覆盖？所占总面积比例是多少？

dev.set(1)

plotPolys(SouthAmBaseProjPS,proj = TRUE,col="wheat1",xlab="longitude",ylab="latitude")

print("Combining Habitat polygons..")

HabL1CB = combinePolys(HabitatL1ProjPS)

HabL2CB = combinePolys(HabitatL2ProjPS)

BrazilCB = combinePolys(BrazilProjPS)

print("...Calculating the UNION of Armadillo Range..")

Armadillo = joinPolys(HabL1CB,HabL2CB,"UNION",maxVert = 1.0e+06)

BrazilArmHab = joinPolys(BrazilCB,Armadillo,"INT",maxVert = 1.0e+06)

BrazilArmHabCB = combinePolys(BrazilArmHab)

attr(BrazilArmHabCB, "projection")

AreaHabInt = calcArea(BrazilArmHabCB,rollup=1)

AllHabitatArea = sum(AreaHabInt[1:length(AreaHabInt[,1]),2])

print(sprintf("Area of Brazil: %g sq km.",BrazilArea))

print(sprintf("Area covered by Armadillo geographic range  (BLUE outline): %g sq km",AllHabitatArea))

addPolys(BrazilProjPS,col="orange")

par(lwd = 2)

addPolys(BrazilArmHabCB,border="blue")

问题五

Armadillo观测点有多少在都在划定的Armadillo分布区之中吗？对分布区内的进行标记，并计算观测点数目

PointsAroundHab = findPolys(SpeciesSightingEventSet,Armadillo)

PointIndexesInsideHab = sort(unique(PointsAroundHab$EID))

PointsInsideHab = SpeciesSightingEventSet[PointIndexesInsideHab,]

allEIDs = seq(length(SpeciesSightingEventSet[,1]))

InsideFlags = !(is.element(allEIDs,PointIndexesInsideHab))

PointIndexesOutsideHab = allEIDs[InsideFlags]

PointsOutsideHab = SpeciesSightingEventSet[PointIndexesOutsideHab,]

print(sprintf("Number of sightings INSIDE geographic range : %d.  Number OUTSIDE geographic range : %d",

length(PointIndexesInsideHab),length(PointIndexesOutsideHab)))

addPoints(PointsInsideHab,col="yellow",pch=19)

addPoints(PointsOutsideHab,col="red",pch=19)

print("Press  to end the demonstration..")

最后修改于 2020-11-09 11:17

阅读(?)评论(0)