v='NLCD2006ImpervLakes_20130212.r'
#Enter Location of NLCD Grid data
Imperv<-'L:/Public/Milstead_Lakes/NLCD2006/nlcd2006_impervious_5-4-11.img'
#Enter loaction of the lake shapefile
MRB1Lakes<-'L:/Public/Milstead_Lakes/WaterbodyDatabase/Shapefiles/MRB1_WBIDLakes.shp'
#load (install) required R packages
if (!'sp' %in% installed.packages()) install.packages('sp');require(sp) #CRS overlay
if (!'maptools' %in% installed.packages()) install.packages('maptools');require(maptools) #readShapePoly
if (!'rgeos' %in% installed.packages()) install.packages('rgeos');require(rgeos) #gBuffer
if (!'rgdal' %in% installed.packages()) install.packages('rgdal');require(rgdal) #readGDAL
if (!'raster' %in% installed.packages()) install.packages('raster');require(raster) #readGDAL
#Projection info
#ESRI USA_Contiguous_Albers_Equal_Area_Conic_USGS_version
AlbersContigUSGS<-CRS('+proj=aea +x_0=0 +y_0=0 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +lat_0=23 +units=m +datum=NAD83')
#ESRI GCS_North_American_1983
NAD83<-CRS("+proj=longlat +datum=NAD83")
#Get the lake shapefile
MRB1Lakes <- readShapePoly(MRB1Lakes, proj4string=NAD83)
#Reproject Lakes to the same CRS as the Population grid
#which(MRB1Lakes$WB_ID==22302965) #27357 index value for lake champlain
MRB1Lakes<-spTransform(MRB1Lakes[-27357,],AlbersContigUSGS) #lake champlain removed
#Associate lake holes (islands) with the correct polygon. This is an important step.
slot(MRB1Lakes, "polygons") <- lapply(slot(MRB1Lakes, "polygons"), checkPolygonsHoles)
#get the NLCD grid data
Imperv<-raster(Imperv)
#image(Imperv)
#extent(Imperv)
#Imperv is the impervious cover for all lower 48 states. Reduce to bbox(MRB1Lakes) + 10km
B<-bbox(MRB1Lakes)
Extent<-c(B[1,1]-10000,B[1,2]+10000,B[2,1]-10000,B[2,2]+10000)
#Get Impervious data for the Largest Lake Buffer + 3 30x30 pixels
Imperv<-crop(Imperv,Extent)
#save the data
File<-paste('L:/Public/Milstead_Lakes/RData/Data_',v,'da',sep='') #file to store data: v is the version (see top of file)
save(Imperv,MRB1Lakes,file=File) #save the last records.
#load("L:/Public/Milstead_Lakes/RData/Data_NLCD2006ImpervLakes_20130131.rda")
################################
##############################
####function to calculate the impervious cover area and percent in the buffer.
CalcImperv<- function(Buf,Imp){
gc() #release unused memory
Mask<-mask(Imp, Buf) #extract impervious cover pixels in the buffer (all others changed to NA)
a<-table(getValues(Mask),useNA='ifany') #the NA are pixels outside of bbox of buffer. Values of "127" are the real NA
a<-na.exclude(data.frame(Percent=as.numeric(names(a)),a)[,-2]) #replace percents stored as factors with values
a$ImpPix<-a$Percent*a$Freq/100 #processing step-convert frequency distribution to number of impervious pixels
MissingPix<-ifelse(max(a$Percent)==127,a[a$Percent==127,'Freq'],0) #number of missing pixels
TotalPix<-sum(a$Freq) #total number of pixels in buffer
PercentNA<-round(MissingPix/TotalPix,3)*100 #percent NA cells in buffer
BufAreaKm2<-round(sum(a$Freq)*30*30/1000000,3) #total area of buffers (based on pixel number)-includes NA (value=127)
a<-subset(a,a$Percent<=100) #remove pixels with missing values (=127)
BufAreaKm2Adj<-round(sum(a$Freq)*30*30/1000000,3) #this doesn't include Pixels with missing data
ImpervAreaKm2<-round(sum(a$ImpPix)*30*30/1000000,3) #calculate the area of impervious cover
PercentImperv<-round(sum(ImpervAreaKm2)/sum(BufAreaKm2Adj),4)*100 #calculate the percent of impervious cover
data.frame(PercentImperv,ImpervAreaKm2,BufAreaKm2,BufAreaKm2Adj,PercentNA)
}
####function to calculate the number of people and density in the buffer.
####################works
Impervious<- function(WB_ID,PlotYN){ #WB_ID ID of lake, PlotYN enter 'Y' to generate figure or 'N' for text only
Start<-Sys.time() #Record Start Time
gc() #release unused memory
Lake<-MRB1Lakes[match(WB_ID,MRB1Lakes$WB_ID),] #select lake
radius<-round(sqrt(Lake$AlbersArea/pi)) #Calculate the approximate lake radius as sqrt(Area/pi)
#Get the buffers
Buf300<-gDifference(gBuffer(Lake,width=300),Lake)
Buf1000<-gDifference(gBuffer(Lake,width=1000),Lake)
Buf2500<-gDifference(gBuffer(Lake,width=2500),Lake)
BufRadius<-gDifference(gBuffer(Lake,width=radius),Lake)
#Select offset and extent for maximum grid if(radius<=2500)B<-bbox(Buf2500) else B<-bbox(BufRadius)
if(radius<=2500)B<-bbox(Buf2500) else B<-bbox(BufRadius)
Extent<-c(B[1,1]-90,B[1,2]+90,B[2,1]-90,B[2,2]+90)
#Get Impervious data for the Largest Lake Buffer + 3 30x30 pixels
Imp<-crop(Imperv,Extent)
#plot grid and buffers ###optional PlotYN='Y' to plot; PlotYN='N' to return text only.
if(PlotYN=='Y'){
image(Imp,main=paste('WB_ID = ',WB_ID))
plot(Buf300,add=T)
plot(Buf1000,add=T)
plot(Buf2500,add=T)
plot(BufRadius,add=T)
title(sub=v)
}
#calculate imperv area and percent for each buffer
a<-data.frame(matrix(nrow=4,ncol=7))
a[1,]<-c(WB_ID,300,CalcImperv(Buf300,Imp))
a[2,]<-c(WB_ID,1000,CalcImperv(Buf1000,Imp))
a[3,]<-c(WB_ID,2500,CalcImperv(Buf2500,Imp))
a[4,]<-c(WB_ID,radius,CalcImperv(BufRadius,Imp))
names(a)<-c('WB_ID','BufWidthM','PercentImperv','ImpervAreaKm2','BufferAreaKm2','BufferAreaKm2Adj','PercentNA')
print(paste(WB_ID,' Time Elapsed = ',round(Sys.time()-Start),' seconds')) #how long did the process take?
print('') #extra line
return(a)
}
############################
#examples
#
Impervious(sample(MRB1Lakes$WB_ID,1),'Y') #random lake without plot of grid and buffers
#Impervious(sample(MRB1Lakes$WB_ID,1),'N') #random lake with results only
#Impervious(4290651,'Y') #big lake
#Impervious(5195222,'Y') #lake on the grid border
#Impervious(1720193,'Y') #lake cannot be buffered with R64 but works in R32. Why?
#Impervious(9326606,'Y') #lake cannot be buffered with R64 but works in R32. Why?
#Impervious(6114592,'Y') #this one looks like a cross section through the small intestines.
#Impervious(Lakes[749,],'Y')
#i<-23;Impervious(MRB1Lakes$WB_ID[i],'Y')
##loop to run a whole bunch of lakes
LakeImperv<-data.frame(matrix(NA,nrow=N*4,ncol=7)) #build data.frame to store the results
names(LakeImperv)<-c('WB_ID','BufWidthM','PercentImperv','ImpervAreaKm2','BufferAreaKm2','BufferAreaKm2Adj','PercentNA')
File<-paste('L:/Public/Milstead_Lakes/RData/',v,'da',sep='') #file to store data: v is the version (see top of file)
#start loop
B<-1 #row number to start processing lakes; this will be 1 to start.
N<-nrow(MRB1Lakes) #last row to process;
#N<-5 #last row to process;
S<-100 #save the work every "S" lakes
Counter<-0 #counter; should be zero to start
for(i in c(B:N)){
LakeImperv[(4*i-3):(i*4),]<-Impervious(MRB1Lakes$WB_ID[i],'N')
Counter<-Counter+1
if (Counter==S) save(LakeImperv,file=File)
if(Counter==S) Counter<-0
}
#End Loop
#add field to distinguish between fixed width and variable (radius) width buffers.
#some radius buffers the same as the fixed width.
LakeImperv$BufType<-c('fixed','fixed','fixed','radius')
#table(LakeImperv$BufType)
#table(LakeImperv[LakeImperv$BufWidthM==300,'BufType'])
#save the data
#save(LakeImperv,file="L:/Public/Milstead_Lakes/RData/NLCD2006ImpervLakes_20130212.rda") #save the last records.
#load("L:/Public/Milstead_Lakes/RData/NLCD2006ImpervLakes_20130212.rda")
#write.table(LakeImperv, file='c:/temp/LakeImperv.csv',row.names=F,sep=',')
# Write data to tblMNKA in fj.mdb
require(RODBC) #Package RODBC must be installed
con <- odbcConnectAccess("C:/Bryan/EPA/Data/WaterbodyDatabase/WaterbodyDatabase.mdb")
FJ <- sqlSave(con,dat=LakeImperv,tablename='LakeNLCD2006ImperviousCover',append=F,rownames=F)
close(con)
##Data Definitions n=28,121 (Lake Champlain excluded due to size and complexity)
#WB_ID: Lake identification number
#BufWidthM: (m) width of lake buffer for calculations
#PercentImperv: (%) Percent impervious cover in buffer = ImpervAreaKm2/BufferAreaKm2Adj/100
#ImpervAreaKm2: (Km2) total impervious cover in the buffer
#BufferAreaKm2: (Km2) Area of entire buffer-number of grid cells*30*30/1000000
#BufferAreaKm2Adj (Km2) Area of Non-NA buffer-number of grid cells with data*30*30/1000000
#PercentNA: (%) (number of NA grid cells)/(total number of grid cells)*100
#BufType: "Fixed" = standard buffer width; "Radius" buffer width = radius of a circle with Area = lake area.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.