Part I: Recursive Data Aquisition and Extraction

In [ ]:
library(maptools)
require(RCurl)
#library(utils)
#URL ftp://prism.nacse.org/ daily monthly

baseurl = "https://aqs.epa.gov/aqsweb/airdata/"

for (year in 1980:2017){
    url = gsub(" ", "", paste(baseurl,"/annual_aqi_by_county_",year,".zip"))
    file = gsub(" ", "", paste("./epa/", year, ".zip"))
    print(paste("downloading",url, " to DEST", file))
    download.file(url, gsub(" ", "", paste(getwd(), "/epa/", year, ".zip", sep = "")))
}


#area <- readShapePoly("ne_10m_parks_and_protected_lands_area.shp")
In [ ]:
url = "ftp://ftp.cpc.ncep.noaa.gov/GIS/us_temp/tmax/"
filenames = getURL(url, ftp.use.epsv = FALSE, dirlistonly = TRUE)
filenames <- strsplit(filenames, "\r\n")
filenames = unlist(filenames)

for (filename in filenames) {
print ("downloading ")
download.file(paste(url, filename, sep = ""), gsub(" ", "", paste(getwd(), "/noaatmax/", filename, sep = "")))
}
In [ ]:
library(curl)
url = "ftp://ftp.cpc.ncep.noaa.gov/GIS/us_temp/tmax/"
h = new_handle(dirlistonly=TRUE)
con = curl(url, "r", h)
tbl = read.table(con, stringsAsFactors=TRUE, fill=TRUE)
close(con)
head(tbl)
for (filename in tbl) {
print ("downloading ")
download.file(paste(url, filename, sep = ""), gsub(" ", "", paste(getwd(), "/noaatmax/", filename, sep = "")))
}
In [ ]:
getwd()
In [ ]:
zipF <- paste(getwd(), "/tmean/", sep = "")
outDir <- paste(getwd(), "/tmean/", sep = "")

for (year in 1895:2016){
unzip(paste(zipF, year, ".zip", sep = ""), exdir=paste(outDir, year, sep=""))
}

Part II: Data Examination and Preparation

All files extracted, now I have to figure out how to deal with the data

First, I should create a single raster for each year

In [ ]:
library(raster)
for (year in 1895:2016){
    testFile = paste(getwd(), "/tmean/", year, "/PRISM_tmean_stable_4kmM2_", year, "_bil.bil", sep="") 
    input_name= raster(testFile)
    output_name= paste(getwd(), "/tmean/", year, "/", year, ".tif", sep="")
    writeRaster(input_name, output_name,format="GTiff",datatype='INT1U',overwrite=TRUE)
    }
In [ ]:
library(raster)
for (year in 1895:2016){
    testFile = paste(getwd(), "/tdmean/", year, "/PRISM_tdmean_stable_4kmM1_", year, "_bil.bil", sep="") 
    input_name= raster(testFile)
    output_name= paste(getwd(), "/tdmean/", year, "/", year, ".tif", sep="")
    writeRaster(input_name, output_name,format="GTiff",datatype='INT1U',overwrite=TRUE)
    }
In [ ]:
library(raster)
for (year in 1895:1980){
    testFile = paste(getwd(), "/precip/", year, "/PRISM_ppt_stable_4kmM2_", year, "_bil.bil", sep="") 
    input_name= raster(testFile)
    output_name= paste(getwd(), "/precip/", year, "/", year, ".tif", sep="")
    writeRaster(input_name, output_name,format="GTiff",datatype='INT1U',overwrite=TRUE)
    }

#switch to M3

for (year in 1981:2016){
    testFile = paste(getwd(), "/precip/", year, "/PRISM_ppt_stable_4kmM3_", year, "_bil.bil", sep="") 
    input_name= raster(testFile)
    output_name= paste(getwd(), "/precip/", year, "/", year, ".tif", sep="")
    writeRaster(input_name, output_name,format="GTiff",datatype='INT1U',overwrite=TRUE)
    }

It would make more sense to make a single RasterBrick have one layer for each year for each variable

In [33]:
library(raster)
#raz = raster(paste(getwd(), "/tmean/2003/2003/2003.tif", sep=""))
#raz
#nrow(raz)
#ncol(raz)
#matr = as.matrix(raz)
#matr[0,0]
xtm <- stack()
for (year in 1895:2016){
    print(paste("adding ", year))
    file = paste(getwd(), "/tmean/", year, "/", year, ".tif", sep="")
    raz = raster(file)
    xtm <- stack( xtm , raz )
}
xtm
#tmean first

[1] "adding  1895"
[1] "adding  1896"
[1] "adding  1897"
[1] "adding  1898"
[1] "adding  1899"
[1] "adding  1900"
[1] "adding  1901"
[1] "adding  1902"
[1] "adding  1903"
[1] "adding  1904"
[1] "adding  1905"
[1] "adding  1906"
[1] "adding  1907"
[1] "adding  1908"
[1] "adding  1909"
[1] "adding  1910"
[1] "adding  1911"
[1] "adding  1912"
[1] "adding  1913"
[1] "adding  1914"
[1] "adding  1915"
[1] "adding  1916"
[1] "adding  1917"
[1] "adding  1918"
[1] "adding  1919"
[1] "adding  1920"
[1] "adding  1921"
[1] "adding  1922"
[1] "adding  1923"
[1] "adding  1924"
[1] "adding  1925"
[1] "adding  1926"
[1] "adding  1927"
[1] "adding  1928"
[1] "adding  1929"
[1] "adding  1930"
[1] "adding  1931"
[1] "adding  1932"
[1] "adding  1933"
[1] "adding  1934"
[1] "adding  1935"
[1] "adding  1936"
[1] "adding  1937"
[1] "adding  1938"
[1] "adding  1939"
[1] "adding  1940"
[1] "adding  1941"
[1] "adding  1942"
[1] "adding  1943"
[1] "adding  1944"
[1] "adding  1945"
[1] "adding  1946"
[1] "adding  1947"
[1] "adding  1948"
[1] "adding  1949"
[1] "adding  1950"
[1] "adding  1951"
[1] "adding  1952"
[1] "adding  1953"
[1] "adding  1954"
[1] "adding  1955"
[1] "adding  1956"
[1] "adding  1957"
[1] "adding  1958"
[1] "adding  1959"
[1] "adding  1960"
[1] "adding  1961"
[1] "adding  1962"
[1] "adding  1963"
[1] "adding  1964"
[1] "adding  1965"
[1] "adding  1966"
[1] "adding  1967"
[1] "adding  1968"
[1] "adding  1969"
[1] "adding  1970"
[1] "adding  1971"
[1] "adding  1972"
[1] "adding  1973"
[1] "adding  1974"
[1] "adding  1975"
[1] "adding  1976"
[1] "adding  1977"
[1] "adding  1978"
[1] "adding  1979"
[1] "adding  1980"
[1] "adding  1981"
[1] "adding  1982"
[1] "adding  1983"
[1] "adding  1984"
[1] "adding  1985"
[1] "adding  1986"
[1] "adding  1987"
[1] "adding  1988"
[1] "adding  1989"
[1] "adding  1990"
[1] "adding  1991"
[1] "adding  1992"
[1] "adding  1993"
[1] "adding  1994"
[1] "adding  1995"
[1] "adding  1996"
[1] "adding  1997"
[1] "adding  1998"
[1] "adding  1999"
[1] "adding  2000"
[1] "adding  2001"
[1] "adding  2002"
[1] "adding  2003"
[1] "adding  2004"
[1] "adding  2005"
[1] "adding  2006"
[1] "adding  2007"
[1] "adding  2008"
[1] "adding  2009"
[1] "adding  2010"
[1] "adding  2011"
[1] "adding  2012"
[1] "adding  2013"
[1] "adding  2014"
[1] "adding  2015"
[1] "adding  2016"

class       : RasterStack 
dimensions  : 621, 1405, 872505, 122  (nrow, ncol, ncell, nlayers)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -125.0208, -66.47917, 24.0625, 49.9375  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
names       : X1895, X1896, X1897, X1898, X1899, X1900, X1901, X1902, X1903, X1904, X1905, X1906, X1907, X1908, X1909, ... 
min values  :     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0,     0, ... 
max values  :    25,    25,    26,    25,    26,    25,    24,    25,    25,    25,    25,    25,    26,    26,    26, ...
In [ ]:
outtm = paste(getwd(), "/tmean/RasterBrick.tif" , sep="")
writeRaster(xtm, filename=outtm, options="INTERLEAVE=BAND", overwrite=TRUE)
In [ ]:
arr = array(xtm[[1:122]])
length(arr)
system.time(arr[500*600])
arr[500*600]
In [ ]:
#now tdmean
xtd <- stack()
for (year in 1895:2016){
    print(paste("adding ", year))
    file = paste(getwd(), "/tdmean/", year, "/", year, ".tif", sep="")
    raz = raster(file)
    xtd <- stack( xtd , raz )
}
xtd
#save file
outtd = paste(getwd(), "/tdmean/RasterBrick.tif" , sep="")
writeRaster(xtd, filename=outtd, overwrite=TRUE)
In [ ]:
#now precip
xpr <- stack()
for (year in 1895:2016){
    print(paste("adding ", year))
    file = paste(getwd(), "/precip/", year, "/", year, ".tif", sep="")
    raz = raster(file)
    xpr <- stack( xpr , raz )
}
xpr
#save file
#outpr = paste(getwd(), "/precip/RasterBrick.tif" , sep="")
#writeRaster(xpr, filename=outpr, overwrite=TRUE)

Now that I have the 3 variables saved into individual TIF rasterbricks, I need to do the same with the 1980-2016 EPA pollution data and make sure the projection is the same. (note: ended up not implementing this in the final design)

In [ ]:
library(raster)
for (year in 1980:2016){
    testFile = paste(getwd(), "/epa/", year, "/annual_aqi_by_county_", year, ".csv", sep="") 
    input_name= raster(testFile)
    output_name= paste(getwd(), "/epa/", year, "/", year, ".tif", sep="")
    writeRaster(input_name, output_name,format="GTiff",datatype='INT1U',overwrite=TRUE)
    }
In [ ]:
library(rgdal)
library(raster)
year = 1990
rawobj <- read.csv(paste(getwd(), "/epa/", year, "/annual_aqi_by_county_", year, ".csv", sep=""), stringsAsFactors = FALSE)

# look at the data structure
rawobj
nrow(rawobj)
In [ ]:
options(geonamesUsername="tbraun96")
library(geonames)
res <- GNsearch(name_equals="Calhoun", q="Alabama", country="US", maxRows=1)
res
In [ ]:
year=1980
options(geonamesUsername="tbraun96")
library(geonames)
colX <- matrix(ncol=1,nrow=867)
colY <- matrix(ncol=1,nrow=867)
rawobj <- read.csv(paste(getwd(), "/epa/", year, "/annual_aqi_by_county_", year, ".csv", sep=""), stringsAsFactors = FALSE)
colX[1:867] = "hello"
colX[3:3] = "yep"
rawobj[25:25,1:1]
for (idx in 1:867){
    res <- GNsearch(name_equals=rawobj[idx:idx, 2:2], q=rawobj[idx:idx,1:1], country="US", maxRows=1)
    if (length(res) != 0){
    colX[idx:idx] = res$lat
    colY[idx:idx] = res$lng
    print(paste("For ", rawobj[idx:idx, 2:2], ",", rawobj[idx:idx,1:1], " we have lat/long of", res$lat, ",", res$lng, sep=""))
        }
    else {
         print(paste("SKIPPING ", rawobj[idx:idx, 2:2], ",", rawobj[idx:idx,1:1]))
         colX[idx:idx] = 0
         colY[idx:idx] = 0
    }
}
tworow <- cbind(colX[1:nrow(rawobj)],colY[1:nrow(rawobj)])
rawobj = cbind(rawobj, tworow)
rawobj
In [ ]:
tworow <- cbind(colX[1:nrow(rawobj)],colY[1:nrow(rawobj)])
rawobj = cbind(rawobj, tworow)
rawobj
In [ ]:
rawobj
In [ ]:
year = 2000
rawobj2 <- read.csv(paste(getwd(), "/epa/", year, "/annual_aqi_by_county_", year, ".csv", sep=""), stringsAsFactors = FALSE)

# look at the data structure
rawobj2

Now I have to convert each table into a raster, then make a raster stack keeping coordinate system homogenous

In [ ]:
library(rgdal)
library(raster)


rbrickref = raster(paste(getwd(), "/precip/RasterBrick.tif" , sep=""))
crs(rbrickref)
colnames(rawobj) <- NULL
rawobj
outf <- SpatialPointsDataFrame(as.matrix(rawobj[,22:23]), rawobj)
outf <- crs("+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs")
In [ ]:

Part III: Statistical Analysis of independent sets

Here, I must take tmean, tdmean, and precip RasterBricks and find a regression for each raster cell by comparing the layers. Thereafter, create an output single-layer raster that extrapolates from the linear regression to 2050 and 2100.

Create Extrapolation for tmean

In [27]:
#Optimized index-value retriever fxn
arrLenPer = 872505 #621x1405
getValues <- function(arr, size, xCoord, yCoord){
        pts = matrix(nrow=1,ncol=size)
        for(year in 1:size){
            pts[year] = arr[((year - 1)*arrLenPer) + (((xCoord-1)*621) + yCoord)]
            pts[year] = replace(pts[year], is.na(pts[year]), "0")
        }
    return(c(pts))
}
In [ ]:
outfile <- writeRaster(xtm, filename='xtm.tif', format="GTiff", overwrite=TRUE,options=c("COMPRESS=LZW"))
print("done saving")
In [138]:
#PrecipBrick <- raster(paste(getwd(), "/precip/RasterBrick.tif" , sep=""))
#TdmeanBrick <- raster(paste(getwd(), "/tdmean/RasterBrick.tif" , sep=""))
#TmeanBrick = xtm
TmeanBrick = array(xtm)
print("Entering loop()")
year=1990
print(paste("adding ", year))
file = paste(getwd(), "/tmean/", year, "/", year, ".tif", sep="")
reference = raster(file)

#tmeanMatrix gives values versus time for a single cell
tmeanExtrapolation2050 = raster(nrow=621,ncol=1405)
tmeanExtrapolation2100 = raster(nrow=621,ncol=1405)
tmeanExtrapolation2050 <- setExtent(tmeanExtrapolation2050, reference, keepres=FALSE, snap=FALSE)
tmeanExtrapolation2100 <- setExtent(tmeanExtrapolation2100, reference, keepres=FALSE, snap=FALSE)
tmeanExtrapolation2050
tmeanExtrapolation2100

for (x in 1:621){
    for(y in 1:1405){
        data = paste("At [", x, ",", y, "] | Completed: ", ((x)/(621))*100, "%", sep="")
        write(data, file=paste(getwd(), "/progress.txt", sep=""), append=FALSE)
        matt = getValues(TmeanBrick, 122, x, y)
        t = lm(formula = matt ~ c(1:122))
        output2050 = as.numeric(t$coefficients[1]) + (as.numeric(t$coefficients[2])*155)
        output2100 = as.numeric(t$coefficients[1]) + (as.numeric(t$coefficients[2])*205)
        tmeanExtrapolation2050[x:x,y:y] = output2050
        tmeanExtrapolation2100[x:x,y:y] = output2100
        
    }
}
print("done")

extrap2050 = paste(getwd(), "/tmean/Extrapolation2050.tif" , sep="")
writeRaster(tmeanExtrapolation2050, filename=extrap2050, format="GTiff", overwrite=TRUE)

extrap2100 = paste(getwd(), "/tmean/Extrapolation2100.tif" , sep="")
writeRaster(tmeanExtrapolation2100, filename=extrap2100, format="GTiff", overwrite=TRUE)

[1] "Entering loop()"
[1] "adding  1990"

class       : RasterLayer 
dimensions  : 621, 1405, 872505  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -125.0208, -66.47917, 24.0625, 49.9375  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 

class       : RasterLayer 
dimensions  : 621, 1405, 872505  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -125.0208, -66.47917, 24.0625, 49.9375  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 

[1] "done"
In [39]:
year=1990
print(paste("adding ", year))
file = paste(getwd(), "/tmean/", year, "/", year, ".tif", sep="")
reference = raster(file)

normalize <- function (raz0){
    out = raster(raz0)
    max = maxValue(raz0)
    arr = array(raz0)
    crs(out) <- crs(raz0)
    for (x in 1:nrow(raz0)){
        for (y in 1:ncol(raz0)){
            data = paste("At [", x, ",", y, "]", sep="")
            write(data, file=paste(getwd(), "/progress.txt", sep=""), append=FALSE)
            if (raz0[x:x,y:y] > 0){
                out[x:x,y:y] = (arr[x*y]/max)* 255
                } 
            else {
                out[x:x,y:y] = 0
            }
        }
    }
    return(out)
}

print("Entering loop")
tmeanExtrapolation2050
hist(tmeanExtrapolation2050)

normed = normalize(tmeanExtrapolation2050)
normed <- setExtent(normed, reference, keepres=FALSE, snap=FALSE)

extrap2050norm = paste(getwd(), "/tmean/Extrapolation2050_norm.tif" , sep="")
writeRaster(normed, filename=extrap2050norm,format="GTiff",datatype='INT1U', overwrite=TRUE)

[1] "adding  1990"
[1] "Entering loop"

class       : RasterLayer 
dimensions  : 621, 1405, 872505  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -125.0208, -66.47917, 24.0625, 49.9375  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : -39.36648, 59.64086  (min, max)
In [147]:
extrap2100 = paste(getwd(), "/tmean/Extrapolation2100.tif" , sep="")
writeRaster(scale(tmeanExtrapolation2100), filename=extrap2100,format="GTiff", overwrite=TRUE)
hist(raster(extrap2100))
plot(raster(extrap2100))

Warning message in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
“11% of the raster cells were used. 100000 values used.”

Create Extrapolation for precip

In [ ]:
PrecipBrick = array(xpr[[92:122]])

#tmeanMatrix gives values versus time for a single cell
precipExtrapolation2050 = raster(nrow=621,ncol=1405)
precipExtrapolation2100 = raster(nrow=621,ncol=1405)
for (x in 1:621){
    for(y in 1:1405){
        data = paste("At [", x, ",", y, "] | Completed: ", ((x)/(621))*100, "%", sep="")
        write(data, file=paste(getwd(), "/progress.txt", sep=""), append=FALSE)
        matt = getValues(PrecipBrick, 30, x, y)
        t = lm(formula = matt ~ c(1:30))
        output2050 = as.numeric(t$coefficients[1]) + (as.numeric(t$coefficients[1])*64)
        output2100 = as.numeric(t$coefficients[1]) + (as.numeric(t$coefficients[1])*114)
        precipExtrapolation2050[x:x,y:y] = output2050
        precipExtrapolation2100[x:x,y:y] = output2100
    }
}

extrap2050 = paste(getwd(), "/precip/Extrapolation2050.tif" , sep="")
writeRaster(precipExtrapolation2050, filename=extrap2050, options="INTERLEAVE=BAND", overwrite=TRUE)

extrap2100 = paste(getwd(), "/precip/Extrapolation2100.tif" , sep="")
writeRaster(precipExtrapolation2100, filename=extrap2100, options="INTERLEAVE=BAND", overwrite=TRUE)
In [ ]:
system.time(PrecipBrick)

Create Extrapolation for tdmean

In [ ]:
TdmeanBrick = array(xtd[[92:122]])

#tmeanMatrix gives values versus time for a single cell
tdmeanExtrapolation2050 = raster(nrow=621,ncol=1405)
tdmeanExtrapolation2100 = raster(nrow=621,ncol=1405)
for (x in 1:621){
    for(y in 1:1405){
        data = paste("At [", x, ",", y, "] | Completed: ", ((x)/(621))*100, "%", sep="")
        write(data, file=paste(getwd(), "/progress.txt", sep=""), append=FALSE)
        matt = getValues(TdmeanBrick, 30, x, y)
        t = lm(formula = matt ~ c(1:30))
        output2050 = as.numeric(t$coefficients[1]) + (as.numeric(t$coefficients[1])*64)
        output2100 = as.numeric(t$coefficients[1]) + (as.numeric(t$coefficients[1])*114)
        tdmeanExtrapolation2050[x:x,y:y] = output2050
        tdmeanExtrapolation2100[x:x,y:y] = output2100
    }
}

extrap2050 = paste(getwd(), "/tdmean/Extrapolation2050.tif" , sep="")
writeRaster(tdmeanExtrapolation2050, filename=extrap2050, options="INTERLEAVE=BAND", overwrite=TRUE)

extrap2100 = paste(getwd(), "/tdmean/Extrapolation2100.tif" , sep="")
writeRaster(tdmeanExtrapolation2100, filename=extrap2100, options="INTERLEAVE=BAND", overwrite=TRUE)

Part IV: Combined Statistical Analysis and Survivability Index (SI)

Now, I will take the output rasters from part III and make a combined raster for 2050 and 2100, showing the survivability index for each region

In [148]:
library(raster)

year=1990
print(paste("adding ", year))
file = paste(getwd(), "/tmean/", year, "/", year, ".tif", sep="")
reference = raster(file)

#I want to normalize the rasters first
normPrecipExtrapolation2050 = raster(paste(getwd(), "/precip/Extrapolation2050.tif" , sep=""))
normPrecipExtrapolation2100 = raster(paste(getwd(), "/precip/Extrapolation2100.tif" , sep=""))
normPrecipExtrapolation2050 <- setExtent(normPrecipExtrapolation2050, reference, keepres=FALSE, snap=FALSE)
normPrecipExtrapolation2100 <- setExtent(normPrecipExtrapolation2100, reference, keepres=FALSE, snap=FALSE)

normTmeanExtrapolation2050 = raster(paste(getwd(), "/tmean/Extrapolation2050.tif" , sep=""))
normTmeanExtrapolation2100 = raster(paste(getwd(), "/tmean/Extrapolation2100.tif" , sep=""))
normTmeanExtrapolation2050 <- setExtent(normTmeanExtrapolation2050, reference, keepres=FALSE, snap=FALSE)
normTmeanExtrapolation2100 <- setExtent(normTmeanExtrapolation2100, reference, keepres=FALSE, snap=FALSE)

normTdmeanExtrapolation2050 = raster(paste(getwd(), "/tdmean/Extrapolation2050.tif" , sep=""))
normTdmeanExtrapolation2100 = raster(paste(getwd(), "/tdmean/Extrapolation2100.tif" , sep=""))
normTdmeanExtrapolation2050 <- setExtent(normTdmeanExtrapolation2050, reference, keepres=FALSE, snap=FALSE)
normTdmeanExtrapolation2100 <- setExtent(normTdmeanExtrapolation2100, reference, keepres=FALSE, snap=FALSE)

#extrap2050 = paste(getwd(), "/tmean/Extrapolation2050.tif" , sep="")
#writeRaster(scale(tmeanExtrapolation2050) +4, filename=extrap2050, format="GTiff",datatype='INT1U', overwrite=TRUE)
hist(normTmeanExtrapolation2050)

print("Done loading and normalizing rasters")

[1] "adding  1990"

Warning message in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
“11% of the raster cells were used. 100000 values used.”
[1] "Done loading and normalizing rasters"
In [149]:
normPrecipExtrapolation2100
arrNormPrecipExtrapolation2050 = scale(array(normPrecipExtrapolation2050))
arrNormTmeanExtrapolation2050 = scale(array(normTmeanExtrapolation2050))
arrNormTdmeanExtrapolation2050 = scale(array(normTdmeanExtrapolation2050))

arrNormPrecipExtrapolation2100 = scale(array(normPrecipExtrapolation2100))
arrNormTmeanExtrapolation2100 = scale(array(normTmeanExtrapolation2100))
arrNormTdmeanExtrapolation2100 = scale(array(normTdmeanExtrapolation2100))

class       : RasterLayer 
dimensions  : 621, 1405, 872505  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -125.0208, -66.47917, 24.0625, 49.9375  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : /data/braunth/precip/Extrapolation2100.tif 
names       : Extrapolation2100 
values      : -2.4871, 6.928239  (min, max)
In [150]:
weightTMean = .5
weightPrecip = .4
weightTdmean = .1
doing2100 = 1

SustainabilityIndex2050 = raster(nrow=621,ncol=1405)
SustainabilityIndex2100 = raster(nrow=621,ncol=1405)

SustainabilityIndex2050 <- setExtent(SustainabilityIndex2050, reference, keepres=FALSE, snap=FALSE)
SustainabilityIndex2100 <- setExtent(SustainabilityIndex2100, reference, keepres=FALSE, snap=FALSE)

arrNormPrecipExtrapolation2050 = array(normPrecipExtrapolation2050)
arrNormTmeanExtrapolation2050 = array(normTmeanExtrapolation2050)
arrNormTdmeanExtrapolation2050 = array(normTdmeanExtrapolation2050)

arrNormPrecipExtrapolation2100 = array(normPrecipExtrapolation2100)
arrNormTmeanExtrapolation2100 = array(normTmeanExtrapolation2100)
arrNormTdmeanExtrapolation2100 = array(normTdmeanExtrapolation2100)

if (doing2100 == 0){
for (x in 1:621){
    for(y in 1:1405){        
        #Rainfall is good, Heat is bad. Score Total: SI(tdmean,precip,tmean) = .4*precip + .4*(1/tmean) + .2*(1/tdmean)
        index = ((x-1)*621) + y
        tmeanVal = arrNormTmeanExtrapolation2050[index]
        tdmeanVal = arrNormTdmeanExtrapolation2050[index]
        
        if (!is.finite(tmeanVal) && is.finite(tmeanVal)){
                SI2050 = (weightPrecip*arrNormPrecipExtrapolation2050[index]*1000/255) - 0 - (weightTdmean*(tdmeanVal*1000/41))
            } 
        else if (is.finite(tmeanVal) && !is.finite(tmeanVal)){
                SI2050 = (weightPrecip*arrNormPrecipExtrapolation2050[index]*1000/255) - (weightTMean*(tmeanVal*1000/60)) - 0
            }
        else if (!is.finite(tmeanVal) && !is.finite(tmeanVal)){
                SI2050 = (weightPrecip*arrNormPrecipExtrapolation2050[index]*1000/255) - 0 - 0
            }
        else {
            SI2050 = (weightPrecip*arrNormPrecipExtrapolation2050[index]*1000/255) - (weightTMean*(tmeanVal*1000/60)) - (weightTdmean*(tdmeanVal*1000/41))
        }
        #SI2100 = (weightPrecip*normPrecipExtrapolation2100[x:x,y:y]) + (weightTMean*(1/normTmeanExtrapolation2100[x:x,y:y])) + (weightTdmean*(1/normTdmeanExtrapolation2100[x:x, y:y]))
        #SI2050[!is.finite(SI2050)] <- 0
        SustainabilityIndex2050[x:x,y:y] = SI2050
        #SustainabilityIndex2100[x:x,y:y] = SI2100
        data = paste("At [", x, ",", y, "] | Completed: ", ((x)/(621))*100, "%", sep="")
        write(data, file=paste(getwd(), "/progress.txt", sep=""), append=FALSE)
            }
        }
    } 

if (doing2100 == 1){
for (x in 1:621){
    for(y in 1:1405){        
        #Rainfall is good, Heat is bad. Score Total: SI(tdmean,precip,tmean) = .4*precip + .4*(1/tmean) + .2*(1/tdmean)
        index = ((x-1)*621) + y
        tmeanVal = arrNormTmeanExtrapolation2100[index]
        tdmeanVal = arrNormTdmeanExtrapolation2100[index]
        
        if (!is.finite(tmeanVal) && is.finite(tmeanVal)){
                SI2100 = (weightPrecip*arrNormPrecipExtrapolation2100[index]*1000/255) - 0 - (weightTdmean*(tdmeanVal*1000/41))
            } 
        else if (is.finite(tmeanVal) && !is.finite(tmeanVal)){
                SI2100 = (weightPrecip*arrNormPrecipExtrapolation2100[index]*1000/255) - (weightTMean*(tmeanVal*1000/60)) - 0
            }
        else if (!is.finite(tmeanVal) && !is.finite(tmeanVal)){
                SI2100 = (weightPrecip*arrNormPrecipExtrapolation2100[index]*1000/255) - 0 - 0
            }
        else {
            SI2100 = (weightPrecip*arrNormPrecipExtrapolation2100[index]*1000/255) - (weightTMean*(tmeanVal*1000/60)) - (weightTdmean*(tdmeanVal*1000/41))
        }
        #SI2050[!is.finite(SI2050)] <- 0
        #SustainabilityIndex2050[x:x,y:y] = SI2050
        SustainabilityIndex2100[x:x,y:y] = SI2100
        data = paste("At [", x, ",", y, "] | Completed: ", ((x)/(621))*100, "%", sep="")
        write(data, file=paste(getwd(), "/progress.txt", sep=""), append=FALSE)
    }
}
    }
In [151]:
extrap2100 = paste(getwd(), "/US_Sustainability_SurvivalIndex_2100.tif", sep="")
writeRaster(SustainabilityIndex2100, filename=extrap2100, format="GTiff", overwrite=TRUE)
#v = raster(paste(getwd(), "/US_Sustainability_SurvivalIndex_2100.tif", sep=""))
#hist(SustainabilityIndex2050)
#hist(scale(SustainabilityIndex2050))
hist(SustainabilityIndex2100)

Finally, output rasters for the sustainability indexes

In [ ]:

Now, figure out how to display nicely

In [31]:
convertToNA <- function (raz){
    out = raster(raz)
    for (x in 1:nrow(raz)){
        for(y in 1:ncol(raz)){
            data = paste("At [", x, ",", y, "] | Completed: ", ((x)/(621))*100, "%", sep="")
        write(data, file=paste(getwd(), "/progress.txt", sep=""), append=FALSE)
            if (raz[x:x,y:y] == 0){
                out[x:x,y:y] = NA
            } else {
                out[x:x, y:y] = raz[x:x,y:y]
            }
        }
    }
    return(out)
}
In [ ]:
partsX = floor(621/20)
partsY = floor(1405/40)
tmp = convertToNA(SustainabilityIndex2050)
tmp
data_matrix <- rasterToPoints(tmp)
head(data_matrix)
obj <- kmeans(data_matrix, 5)
plot(data_matrix, col = obj$cluster)
In [ ]:
#Do the same thing as above except normalize once again

extrap2050 = paste(getwd(), "/US_Sustainability_SurvivalIndex_2050_norm.tif" , sep="")
writeRaster(normImage(SustainabilityIndex2050, norm=TRUE), filename=extrap2050, options="INTERLEAVE=BAND", overwrite=TRUE)

extrap2100 = paste(getwd(), "/US_Sustainability_SurvivalIndex_2100_norm.tif" , sep="")
writeRaster(normImage(SustainabilityIndex2100, norm=TRUE), filename=extrap2100, options="INTERLEAVE=BAND", overwrite=TRUE)