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

[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  :    19,    19,    20,    19,    20,    20,    18,    19,    19,    19,    19,    19,    19,    20,    19, ...
In [47]:
outfile <- writeRaster(xtd, filename='xtd.tif', format="GTiff", overwrite=TRUE,options=c("COMPRESS=LZW"))
print("done saving")

[1] "done saving"
In [33]:
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 [48]:
TdmeanBrick = array(xtd)
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
tdmeanExtrapolation2050 = raster(nrow=621,ncol=1405)
tdmeanExtrapolation2100 = raster(nrow=621,ncol=1405)
tdmeanExtrapolation2050 <- setExtent(tdmeanExtrapolation2050, reference, keepres=FALSE, snap=FALSE)
tdmeanExtrapolation2100 <- setExtent(tdmeanExtrapolation2100, reference, keepres=FALSE, snap=FALSE)
tdmeanExtrapolation2050
tdmeanExtrapolation2100

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(), "/progresstd.txt", sep=""), append=FALSE)
        matt = getValues(TdmeanBrick, 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)
        tdmeanExtrapolation2050[x:x,y:y] = output2050
        tdmeanExtrapolation2100[x:x,y:y] = output2100
    }
}

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

extrap2100 = paste(getwd(), "/tdmean/Extrapolation2100.tif" , sep="")
writeRaster(tdmeanExtrapolation2100, 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 [50]:
extrap2100 = paste(getwd(), "/tdmean/Extrapolation2100.tif" , sep="")
writeRaster(scale(tdmeanExtrapolation2100), 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.”
In [30]:
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(), "/progresstd.txt", sep=""), append=FALSE)
            if (!is.na(raz0[x:x,y:y])){
            if (raz0[x:x,y:y] > 0){
                out[x:x,y:y] = (arr[x*y]/max)* 255
                } 
            else {
                out[x:x,y:y] = 0
            }
                } else {
                out[x:x,y:y] = 0
            }
        }
    }
    return(out)
}

tdmeanExtrapolation2050
hist(tdmeanExtrapolation2050)

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

extrap2050norm = paste(getwd(), "/tdmean/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      : -1.418465, 2.077419  (min, max)
In [17]:
library(raster)
for (year in 2016:2016){
    testFile = paste(getwd(), "/tmean/", year, "/", year, ".tif", sep="") 
    input_name= raster(testFile)
    plot(input_name)
    input_name
    #writeRaster(input_name, output_name,format="GTiff",datatype='INT1U',overwrite=TRUE)
    }
In [24]:
year=1990
print(paste("adding ", year))
file = paste(getwd(), "/tmean/", year, "/", year, ".tif", sep="")
extent(raz)
extent(normed)

normed = setExtent(normed, raz, keepres=FALSE, snap=FALSE)
extent(normed)

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

[1] "adding  1990"

class       : Extent 
xmin        : -125.0208 
xmax        : -66.47917 
ymin        : 24.0625 
ymax        : 49.9375 
class       : Extent 
xmin        : -125.0208 
xmax        : -66.47917 
ymin        : 24.0625 
ymax        : 49.9375 
class       : Extent 
xmin        : -125.0208 
xmax        : -66.47917 
ymin        : 24.0625 
ymax        : 49.9375
In [ ]: