In [11]:
library(raster)
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

[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  :    23,    33,    57,    18,    17,    24,    45,    37,     8,    11,    25,    68,    55,    52,    52, ... 
max values  :   255,   255,   255,   255,   255,   255,   255,   255,   255,   255,   255,   255,   255,   255,   255, ...
In [21]:
outfile <- writeRaster(xpr, filename='xpr.tif', format="GTiff", overwrite=TRUE,options=c("COMPRESS=LZW"))
print("done saving")

[1] "done saving"
In [10]:
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 [22]:
PrecipBrick = array(xpr)
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
precipExtrapolation2050 = raster(nrow=621,ncol=1405)
precipExtrapolation2100 = raster(nrow=621,ncol=1405)
precipExtrapolation2050 <- setExtent(precipExtrapolation2050, reference, keepres=FALSE, snap=FALSE)
precipExtrapolation2100 <- setExtent(precipExtrapolation2100, reference, keepres=FALSE, snap=FALSE)
precipExtrapolation2050
precipExtrapolation2100

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(), "/progressprecip.txt", sep=""), append=FALSE)
        matt = getValues(PrecipBrick, 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)
        precipExtrapolation2050[x:x,y:y] = output2050
        precipExtrapolation2100[x:x,y:y] = output2100
    }
}

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

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

[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
In [25]:
extrap2100 = paste(getwd(), "/precip/Extrapolation2100.tif" , sep="")
writeRaster(scale(precipExtrapolation2100), filename=extrap2100, format="GTiff", overwrite=TRUE)
hist(raster(paste(getwd(), "/precip/Extrapolation2100.tif" , sep="")))
plot(raster(extrap2100))

Warning message in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
“11% of the raster cells were used. 100000 values used.”
In [9]:
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(), "/progressprecip.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)
}

precipExtrapolation2050
hist(precipExtrapolation2050)

normed = normalize(precipExtrapolation2050)
hist(normed)
plot(normed)
normed <- setExtent(normed, reference, keepres=FALSE, snap=FALSE)

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

[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 
data source : in memory
names       : layer 
values      : -341.102, 483.2145  (min, max)
In [ ]: