Critical Infrastructure Vulnerability Analysis

Project Objective:

To determine if critical infrastructure in East Central Florida is vulnerable by finding how much they are clustered around each other. As the final project to my Algorithms in Geospatial Intelligence class, I structured the project and the methods around using the R coding language. 

Research Question: 

How vulnerable is critical infrastructure in East Central Florida?

Background:

Critical Infrastructure – Assets, systems, and networks that are considered so vital to the US that their incapacitation would have a debilitating effect on security, national economic security, national public health or safety, or any combination thereof (The Cybersecurity and Infrastructure Security Agency)

Themes from Past Research

• Building critical infrastructure resilience is important to prevent a failure in the system from cascading into a bigger problem.

• The clustering of critical infrastructure creates vulnerability from natural disasters or man-made attacks

• Florida’s critical infrastructure is endangered by natural disasters like hurricanes, which can disrupt energy and transportation networks via flooding

Importance

• Gives policymaker an evaluation of critical infrastructure’s vulnerability in East Central Florida, paving the way for improvements to be made to enhance resiliency before natural disaster strikes

Hypotheses:

Null Hypothesis (H0) - As a collective, critical infrastructure in East Central Florida is not vulnerable

Alternative Hypothesis (H1) - As a collective, critical infrastructure in East Central Florida is vulnerable

Data

Region of Interest: East Central Florida (Brevard, Lake, Orange, Osceola, Seminole, Volusia)

• County shapefiles were found on the Florida Department of Transportation's Open Data Hub.

Critical Infrastructure points were taken from open sources that covered:

• Bridges

• Hospitals

• Airports

• Power Plants

• Water Treatment Facilities

Methods:

1. Imported my compiled point pattern and shapefile into R and projected them both to NAD 1983 UTM Zone 17N

2. Completed Ghat w/ Monte Carlo simulation envelopes to check nearest-neighbor scale clustering in the point pattern

3. Completed Lhat w/ Monte Carlo simulation envelopes to check clustering at different scales

4. Ran Clark-Evans R as a robustness check since Clark-Evans R details a global summary of overall clustering intensity

Results

Ghat: Observed Ghat line went above the uppermost simulation envelope at short distances, signaling statistically significant clustering in the observed point pattern

Lhat: Observed Lhat line went above the uppermost simulation envelope at all distances, signaling statistically significant clustering in the observed point pattern

Clark-Evans R: The expected mean nearest neighbor distance divided by the observed returns a value less than 1. Since 0.5852 is far from 1, we determine that there is a strong clustering pattern in the observed point  pattern. 

Based on the evidence from the three tests, we reject the null hypothesis. Critical infrastructure in East Central Florida is vulnerable.

Further Research Opportunities:

• Has critical infrastructure in East Central Florida gotten more clustered over time?

• Was failing critical infrastructure during Hurricane Milton clustered in a single area?

R Code:

Importing CSV and Shapefile, projecting both to NAD 1983 UTM Zone 17N, and setting up polygon

> Coords<-read.csv("C:\\...\\Documents\\GEOG665_StudySetProjectedCoords.csv")

> XYPoints<-as.points(Coords)

> XYPoints<-as.data.frame(XYPoints)

> XYPoints<-st_as_sf(XYPoints,coords=c("arx","ary"),crs=26917)

> CI.pts<-st_coordinates(XYPoints)

> ECF_Shapefile<-st_read("E:\\...\\County Boundaries\\East Central Florida\\FloridaCountyBoundarieswithFDOTDistricts.shp")

> ECF_Shapefile<-st_transform(ECF_Shapefile,26917)

> ECF.area<-sum(st_area(ECF_Shapefile))

> ECF.area

16279395655 [m^2]

> ECF_union<-st_union(ECF_Shapefile)

> ECF_sp <- as(ECF_union, "Spatial")

> coords_list <- lapply(ECF_sp@polygons, function(x) x@Polygons[[1]]@coords)

> ECF_polygon <- list(x = coords_list[[1]][,1], y = coords_list[[1]][,2])

> ECF_poly_matrix<-cbind(ECF_polygon$x,ECF_polygon$y)

Ghat

Ghat_Monte_Carlo<-function(n.pts,Polygon,ntrials=99){#gives simulation envelopes

  pts.nndist=nndistG(n.pts)

  max.nndist=max(pts.nndist$dist)

  dist.bin=seq(0,max.nndist,max.nndist/100)

  lambda=npts(n.pts)/areapl(Polygon)

  expected.ghat=1-exp(-lambda*pi*d^2)

  n.ghat=Ghat(n.pts,dist.bin)

  plot(dist.bin,n.ghat,type="l",col='green',lwd=2,xlab="distance",ylab="Ghat")

  lines(d,expected.ghat,col='blue',lwd=2)

  hold.mat=matrix(nrow=length(dist.bin),ncol=ntrials,NA)

  for (i in 1:ntrials) {

    trand=gen(Polygon,npts(n.pts))

    hold.mat[,i]=Ghat(trand,dist.bin)

  }

  trand.max=apply(hold.mat,1,max)

  trand.min=apply(hold.mat,1,min)

  trand.mean=apply(hold.mat,1,mean)

  sim_ghat<-cbind(trand.min,trand.max,trand.max)

  lines(dist.bin,sim_ghat[,1],lty=2,lwd=2)

  lines(dist.bin,sim_ghat[,3],lty=2,lwd=2)

  return(legend("bottomright",legend=c("Observed Ghat","Expected Ghat", "Simulation Envelopes (Max,Min)"),col=c("green","blue","black"),lty=c(1,1,2),lwd=c(2,2,2),bty="n"))

}

> Ghat_Monte_Carlo(CI.pts,ECF_poly_matrix,99)

Lhat

> Lhat_Monte_Carlo=function(n.pts,Polygon,ntrials=99){

+  n.pts<-n.pts/1000

+  Polygon<-Polygon/1000

+  minx=min(n.pts[,1])

+  maxx=max(n.pts[,1])

+  miny=min(n.pts[,2])

+  maxy=max(n.pts[,2])

+  num.points<-npts(n.pts)

+  L<-min(maxx-minx,maxy-miny) 

+  d<-seq(0,L,L/100)

+  ou.khat=khat(n.pts,Polygon,d)

+  lhat.out=sqrt(ou.khat/pi)-d

+  plot(d,lhat.out,type="l",col='green',lwd=2,xlab="Distance",ylab="Lhat",main="Lhat with Monte Carlo Envelopes")

+  lines(d,rep(0,length(d)),col='blue',lwd=2)

+  tlist=Kenv.csr(num.points,Polygon,ntrials,d,quiet=T)

+  tlist$lower=sqrt(tlist$lower/pi)-d

+  tlist$upper=sqrt(tlist$upper/pi)-d

+  lines(d,tlist$lower,lty=2,lwd=2)

+  lines(d,tlist$upper,lty=2,lwd=2)

+  return(legend("bottomright",legend=c("Observed Lhat","Expected Lhat","Simulation Envelopes (Max,Min)"),col=c("green","blue","black"),lty=c(1,1,2),lwd=c(2,2,2),bty="n"))

> Lhat_Monte_Carlo(CI.pts,ECF_poly_matrix,20)

Clark-Evans R

> CI.nndist<-nndistG(CI.pts)

> mean(CI.nndist$dists)

[1] 914.1041

> CI.npts<-npts(CI.pts)

> R.exp=1/(2*sqrt(CI.npts/ECF.area))

> R.exp

1562.036 [m]

> Final_R=914.1041/1562.036

> Final_R

[1] 0.5852004