Workshop 3
Natasha Langton
2 April 2019
library(sp)
library(raster)
library(ncdf4)
library(rgdal)## Warning: package 'rgdal' was built under R version 3.5.3## rgdal: version: 1.4-3, (SVN revision 828)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/Natasha/Documents/R/win-library/3.5/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/Natasha/Documents/R/win-library/3.5/rgdal/proj
## Linking to sp version: 1.3-1library(RColorBrewer)
library(lattice)
library(latticeExtra)
#library(OpenStreetMap)
library(reshape2)
library(maps)clear variables and any open plots before starting
graphics.off()
rm(list=ls())set wd
setwd("~/College/Msc Climate Change/GY667/Workshop 3")path <- file.path(getwd(),"data");read coastline using readOGR as a Spatial Data Frame part of the rgdal package higher resolution coastlines available but not needed for this application
world.coast <- readOGR(dsn = file.path(path,"ne_110m_coastline"), layer = "ne_110m_coastline" )## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\Natasha\Documents\College\Msc Climate Change\GY667\Workshop 3\data\ne_110m_coastline", layer: "ne_110m_coastline"
## with 134 features
## It has 3 fields
## Integer64 fields read as strings: scalerankconvert to the correct projections
world.coast <- spTransform(world.coast, CRS( "+proj=longlat +datum=WGS84" ) )task 1:
#load HadISST_sst.nc
nc <- nc_open(file.path(getwd(), "HADISST_sst.nc"))
lat <- ncvar_get(nc,"latitude")
lon <- ncvar_get(nc,"longitude")
sst <- ncvar_get(nc, "sst")
date <- ncvar_get(nc, "time")Get units, always changeable with time
tunits<-ncatt_get(nc,"time",attname="units")tunits = “days since 1850-1-1 00:00:00”
tustr<-strsplit(tunits$value, " ")Check you are happy with the data time origin
date<-as.character(as.Date(date,origin=unlist(tustr)[3]))replace missing values with NA, many datasets have different fillvalues
fillvalue <- ncatt_get(nc,"sst","_FillValue")
sst[sst==fillvalue$value] <- NA
missvalue <- ncatt_get(nc,"sst","missing_value")
sst[sst==missvalue$value] <- NA
sst[sst==-1000] <- NAwe’re interested in patterns of global variability so we need to remove the global mean at each step but first investigate the global mean
use annual averages with aggregate and using colMeans and rowMeans you could alternatively use apply as in the first workshop
year <- format(as.Date(date, format="%Y-%m-%d"),"%Y")
gmean <- colMeans(sst, na.rm = TRUE, dims=2)
annmean <- aggregate(gmean,by=list(year),FUN=mean,na.rm=TRUE)alternatively comment out the following two lines, what difference does excluding NA values make?
avsst = rowMeans(sst,na.rm=TRUE,dims=2)
head(avsst)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] -1.799712 -1.8 -1.8 -1.799470 -1.794967 -1.792652 -1.791727 -1.792087
## [2,] -1.800000 -1.8 -1.8 -1.799746 -1.796744 -1.789094 -1.791494 -1.791266
## [3,] -1.799768 -1.8 -1.8 -1.799362 -1.796093 -1.793520 -1.791205 -1.790556
## [4,] -1.799773 -1.8 -1.8 -1.799630 -1.797190 -1.792647 -1.788093 -1.791084
## [5,] -1.799785 -1.8 -1.8 -1.799896 -1.795571 -1.787489 -1.793135 -1.790110
## [6,] -1.799801 -1.8 -1.8 -1.799905 -1.799620 -1.793820 -1.791726 -1.791580
## [,9] [,10] [,11] [,12] [,13] [,14] [,15]
## [1,] -1.790149 -1.786665 -1.780915 -1.773490 -1.762206 -1.741153 -1.719786
## [2,] -1.791421 -1.787862 -1.782248 -1.774436 -1.764528 -1.740194 -1.715657
## [3,] -1.790178 -1.783250 -1.780777 -1.775022 -1.763293 -1.730552 -1.700452
## [4,] -1.789861 -1.784812 -1.772676 -1.771306 -1.761183 -1.736301 -1.701224
## [5,] -1.792168 -1.781883 -1.776275 -1.771795 -1.755532 -1.733240 -1.693643
## [6,] -1.792786 -1.780379 -1.776259 -1.771065 -1.752386 -1.722562 -1.694405
## [,16] [,17] [,18] [,19] [,20] [,21]
## [1,] -1.666631 -1.555555 -1.441381 -1.1743229 -0.7754281 -0.6493827
## [2,] -1.672443 -1.572148 -1.379343 -1.0760637 -0.8533030 -0.4694451
## [3,] -1.668149 -1.564931 -1.374210 -1.0966351 -0.7833820 -0.4159841
## [4,] -1.667858 -1.532012 -1.352431 -0.9730221 -0.6486905 -0.2932194
## [5,] -1.655646 -1.536851 -1.335464 -0.8351600 -0.6034487 -0.1584155
## [6,] -1.654446 -1.534960 -1.126895 -0.7888845 -0.3697352 -0.1009892
## [,22] [,23] [,24] [,25] [,26] [,27] [,28]
## [1,] NaN NaN NaN 1.0249377 1.126673 1.262509 1.773419
## [2,] NaN NaN NaN 1.0452419 1.379340 1.618507 1.822994
## [3,] -0.35531499 NaN NaN 0.9135552 1.258294 1.734846 1.906655
## [4,] -0.04176329 NaN NaN NaN 1.302797 1.760991 1.973075
## [5,] -0.08089632 NaN NaN NaN 1.355508 1.804370 2.003524
## [6,] -0.03962538 -0.06106675 NaN NaN 1.620295 1.615157 2.027330
## [,29] [,30] [,31] [,32] [,33] [,34] [,35]
## [1,] 2.644969 3.381521 3.858063 4.336045 4.717195 4.937529 5.121027
## [2,] 2.640299 3.343127 3.856099 4.347613 4.757762 4.992551 5.178441
## [3,] 2.605906 3.285104 3.853700 4.337926 4.789467 5.051107 5.245836
## [4,] 2.580001 3.209619 3.823459 4.288603 4.799861 5.101757 5.308497
## [5,] 2.465389 3.125613 3.760868 4.249798 4.791644 5.134830 5.358847
## [6,] 2.367835 3.046399 3.676972 4.190188 4.757562 5.136606 5.383885
## [,36] [,37] [,38] [,39] [,40] [,41] [,42]
## [1,] 5.228368 5.276693 5.383786 5.588450 5.834319 6.089891 6.392899
## [2,] 5.282075 5.332977 5.448545 5.654996 5.893275 6.142618 6.446611
## [3,] 5.352025 5.406577 5.527186 5.732940 5.963246 6.205951 6.509534
## [4,] 5.415242 5.471367 5.595007 5.799395 6.024938 6.265406 6.569547
## [5,] 5.466341 5.523685 5.649658 5.852417 6.075754 6.317786 6.623745
## [6,] 5.499960 5.561589 5.691114 5.894467 6.116872 6.360271 6.668809
## [,43] [,44] [,45] [,46] [,47] [,48] [,49]
## [1,] 6.807126 7.395176 8.198454 9.182922 10.28493 11.46250 12.66317
## [2,] 6.863599 7.448278 8.239780 9.208206 10.29554 11.46731 12.67373
## [3,] 6.930888 7.519525 8.308837 9.270089 10.34655 11.50894 12.71464
## [4,] 6.994052 7.589315 8.382902 9.340843 10.40830 11.55934 12.75655
## [5,] 7.051349 7.656564 8.459444 9.415292 10.47420 11.61128 12.79043
## [6,] 7.099608 7.711499 8.520291 9.473809 10.52468 11.64938 12.81133
## [,50] [,51] [,52] [,53] [,54] [,55] [,56]
## [1,] 13.82411 14.92240 15.94492 16.87619 17.71784 18.47638 19.17034
## [2,] 13.84323 14.94294 15.95125 16.85455 17.67240 18.42155 19.11552
## [3,] 13.88537 14.97584 15.96337 16.84375 17.64730 18.39475 19.09190
## [4,] 13.92118 14.99712 15.96288 16.82730 17.62625 18.37675 19.07501
## [5,] 13.94158 15.00250 15.95163 16.80919 17.61009 18.36252 19.05641
## [6,] 13.94616 14.99220 15.93050 16.78829 17.59514 18.34791 19.03604
## [,57] [,58] [,59] [,60] [,61] [,62] [,63]
## [1,] 19.83273 20.48815 21.13088 21.77575 22.45981 23.17053 23.84759
## [2,] 19.78132 20.44299 21.09587 21.74882 22.43300 23.14096 23.81549
## [3,] 19.75929 20.42428 21.08470 21.74278 22.42178 23.12099 23.79121
## [4,] 19.73970 20.40654 21.07386 21.73468 22.40468 23.09026 23.75846
## [5,] 19.71551 20.38603 21.06221 21.72662 22.38800 23.05630 23.71929
## [6,] 19.69304 20.37081 21.05814 21.72872 22.38244 23.03004 23.67799
## [,64] [,65] [,66] [,67] [,68] [,69] [,70]
## [1,] 24.44503 24.94911 25.36650 25.72611 26.04796 26.32445 26.55173
## [2,] 24.40719 24.90331 25.31559 25.67516 25.99783 26.27334 26.50220
## [3,] 24.38067 24.87531 25.28527 25.64079 25.95887 26.23155 26.46226
## [4,] 24.35744 24.85826 25.26346 25.60874 25.92117 26.19730 26.43565
## [5,] 24.32973 24.83760 25.23590 25.56941 25.87904 26.16596 26.41685
## [6,] 24.29044 24.80547 25.20258 25.52846 25.83635 26.13248 26.39359
## [,71] [,72] [,73] [,74] [,75] [,76] [,77]
## [1,] 26.75948 26.97560 27.17508 27.31499 27.39569 27.44509 27.50066
## [2,] 26.71651 26.93879 27.13534 27.26364 27.33830 27.39580 27.46339
## [3,] 26.68056 26.90190 27.08964 27.20817 27.28132 27.34564 27.42130
## [4,] 26.65797 26.87157 27.04160 27.14863 27.22409 27.29603 27.37678
## [5,] 26.64216 26.84217 26.98924 27.08618 27.16728 27.24643 27.33084
## [6,] 26.61975 26.80794 26.93912 27.03194 27.12083 27.20686 27.29251
## [,78] [,79] [,80] [,81] [,82] [,83] [,84]
## [1,] 27.58782 27.70378 27.83576 27.97417 28.12420 28.27210 28.39186
## [2,] 27.55613 27.67277 27.80587 27.94860 28.10446 28.25597 28.37095
## [3,] 27.52063 27.64253 27.77900 27.92458 28.08102 28.22859 28.33871
## [4,] 27.48196 27.61461 27.76514 27.92178 28.07634 28.21087 28.30764
## [5,] 27.44234 27.58792 27.75747 27.93150 28.08726 28.20787 28.28541
## [6,] 27.40798 27.56726 27.75564 27.94323 28.10084 28.21430 28.27829
## [,85] [,86] [,87] [,88] [,89] [,90] [,91]
## [1,] 28.47415 28.51587 28.52326 28.51095 28.49274 28.48855 28.52744
## [2,] 28.43924 28.46788 28.47116 28.46078 28.44407 28.44197 28.48622
## [3,] 28.40339 28.42512 28.41499 28.38883 28.36379 28.36995 28.43683
## [4,] 28.36446 28.37903 28.35281 28.30614 28.27174 28.29055 28.38962
## [5,] 28.32670 28.33317 28.29583 28.23502 28.19494 28.22380 28.35119
## [6,] 28.30480 28.30233 28.25530 28.18192 28.13556 28.16775 28.30796
## [,92] [,93] [,94] [,95] [,96] [,97] [,98]
## [1,] 28.61702 28.74300 28.88531 29.02458 29.15066 29.24250 29.27339
## [2,] 28.58222 28.71505 28.86327 28.99984 29.11747 29.20710 29.24266
## [3,] 28.55614 28.70197 28.85207 28.98083 29.08622 29.16990 29.20980
## [4,] 28.54268 28.70200 28.84349 28.95996 29.06033 29.14489 29.19069
## [5,] 28.53877 28.71141 28.83930 28.94378 29.04706 29.13829 29.19097
## [6,] 28.51316 28.69661 28.82501 28.92744 29.03195 29.12959 29.19383
## [,99] [,100] [,101] [,102] [,103] [,104] [,105]
## [1,] 29.24505 29.18890 29.11939 29.01377 28.85041 28.61887 28.31397
## [2,] 29.22484 29.18518 29.13362 29.04275 28.88339 28.64399 28.33023
## [3,] 29.20218 29.17470 29.13467 29.05130 28.89371 28.65287 28.34189
## [4,] 29.19297 29.17226 29.13106 29.04408 28.88879 28.66180 28.37235
## [5,] 29.20288 29.18434 29.13087 29.03048 28.87656 28.67322 28.41692
## [6,] 29.21767 29.19998 29.13420 29.02181 28.86612 28.67588 28.44168
## [,106] [,107] [,108] [,109] [,110] [,111] [,112]
## [1,] 27.95036 27.55969 27.17528 26.79081 26.35562 25.85138 25.31002
## [2,] 27.96680 27.58841 27.22257 26.84872 26.40623 25.88464 25.32946
## [3,] 27.98724 27.61537 27.24897 26.86978 26.41777 25.88461 25.32405
## [4,] 28.03292 27.66220 27.28093 26.87706 26.40647 25.86653 25.31046
## [5,] 28.09580 27.72420 27.31910 26.87248 26.37227 25.82916 25.28718
## [6,] 28.13407 27.76081 27.33769 26.86188 26.34389 25.80253 25.26955
## [,113] [,114] [,115] [,116] [,117] [,118] [,119]
## [1,] 24.77544 24.26276 23.74597 23.21039 22.67361 22.14416 21.59859
## [2,] 24.79243 24.28648 23.77252 23.23126 22.68172 22.13107 21.55917
## [3,] 24.78959 24.28694 23.77442 23.22981 22.66252 22.08592 21.49868
## [4,] 24.78290 24.27734 23.75738 23.20474 22.62644 22.04467 21.46778
## [5,] 24.77164 24.26021 23.72834 23.16610 22.58412 22.01433 21.46285
## [6,] 24.75463 24.23170 23.69317 23.13501 22.56103 22.00482 21.47340
## [,120] [,121] [,122] [,123] [,124] [,125] [,126]
## [1,] 21.00727 20.39410 19.82654 19.31977 18.84329 18.38774 17.95152
## [2,] 20.94474 20.32369 19.76499 19.26545 18.77950 18.30185 17.84551
## [3,] 20.88716 20.28231 19.73482 19.22532 18.71058 18.19802 17.71666
## [4,] 20.88108 20.30038 19.75273 19.21278 18.65333 18.10097 17.59691
## [5,] 20.90624 20.34335 19.78324 19.20708 18.60652 18.02028 17.49466
## [6,] 20.93432 20.37059 19.78683 19.18060 18.55815 17.95510 17.41028
## [,127] [,128] [,129] [,130] [,131] [,132] [,133]
## [1,] 17.51797 17.06518 16.60105 16.15039 15.67637 15.07315 14.25427
## [2,] 17.40998 16.98034 16.54860 16.11970 15.66300 15.08177 14.29197
## [3,] 17.27593 16.85351 16.43094 16.01410 15.57471 15.02137 14.27896
## [4,] 17.14795 16.72215 16.29981 15.89102 15.46969 14.94928 14.26540
## [5,] 17.02976 16.59160 16.16698 15.76438 15.35857 14.87204 14.25147
## [6,] 16.92409 16.47031 16.04095 15.64127 15.24576 14.78554 14.21496
## [,134] [,135] [,136] [,137] [,138] [,139] [,140]
## [1,] 13.27339 12.30132 11.46592 10.75858 10.08537 9.416843 8.891008
## [2,] 13.34771 12.40458 11.58651 10.89705 10.24591 9.596250 9.073308
## [3,] 13.39439 12.50756 11.73678 11.08742 10.47371 9.851744 9.327616
## [4,] 13.45987 12.64642 11.93071 11.32665 10.75543 10.165143 9.634758
## [5,] 13.53148 12.79128 12.12069 11.54936 11.01216 10.457214 9.929723
## [6,] 13.56059 12.87762 12.23904 11.68490 11.17282 10.651357 10.134795
## [,141] [,142] [,143] [,144] [,145] [,146] [,147]
## [1,] 8.423951 8.030039 7.617577 7.159713 6.653107 6.096617 5.477840
## [2,] 8.596709 8.180797 7.740390 7.254704 6.721100 6.135208 5.488662
## [3,] 8.824090 8.361171 7.873391 7.350591 6.787400 6.176583 5.504535
## [4,] 9.085650 8.554732 8.009507 7.441781 6.839059 6.197486 5.503200
## [5,] 9.330962 8.727171 8.125147 7.512736 6.867696 6.189107 5.474352
## [6,] 9.509652 8.851994 8.201692 7.549889 6.870418 6.158705 5.422975
## [,148] [,149] [,150] [,151] [,152] [,153] [,154]
## [1,] 4.792800 4.051160 3.288124 2.563124 1.893136 1.308794 0.7784685
## [2,] 4.782238 4.025948 3.254650 2.524728 1.874019 1.306799 0.7908744
## [3,] 4.772963 3.998688 3.217939 2.484215 1.842973 1.285715 0.7773835
## [4,] 4.754645 3.970628 3.186019 2.449848 1.811585 1.254575 0.7492272
## [5,] 4.721331 3.938311 3.154519 2.412405 1.770459 1.210653 0.7045036
## [6,] 4.665592 3.886683 3.107264 2.364477 1.726514 1.158537 0.6669740
## [,155] [,156] [,157] [,158] [,159] [,160]
## [1,] 0.2631246 -0.2464980 -0.6583665 -0.9745090 -1.194044 -1.323445
## [2,] 0.3058383 -0.2135148 -0.6278392 -0.9716997 -1.188790 -1.294794
## [3,] 0.3033367 -0.2135408 -0.6210952 -0.9769241 -1.193549 -1.297581
## [4,] 0.2851924 -0.2158019 -0.6273059 -0.9807599 -1.199231 -1.302827
## [5,] 0.2491880 -0.2386569 -0.6398985 -0.9791211 -1.205520 -1.308180
## [6,] 0.2311133 -0.2632437 -0.6598558 -0.9766402 -1.209952 -1.313924
## [,161] [,162] [,163] [,164] [,165] [,166] [,167]
## [1,] -1.360775 -1.376471 -1.369580 -1.314345 -1.287331 -1.286518 -1.291519
## [2,] -1.361393 -1.379028 -1.358846 -1.315172 -1.293562 -1.298143 -1.314858
## [3,] -1.362431 -1.383943 -1.363874 -1.323093 -1.301877 -1.301635 -1.323908
## [4,] -1.365636 -1.388522 -1.377493 -1.348421 -1.313195 -1.307584 -1.336111
## [5,] -1.377207 -1.395174 -1.399367 -1.377362 -1.341177 -1.325124 -1.349523
## [6,] -1.380703 -1.407039 -1.420193 -1.411699 -1.372397 -1.343477 -1.360978
## [,168] [,169] [,170] [,171] [,172] [,173] [,174] [,175] [,176]
## [1,] -1.232284 NaN NaN NaN NaN NaN NaN NaN NaN
## [2,] -1.264675 NaN NaN NaN NaN NaN NaN NaN NaN
## [3,] -1.289101 NaN NaN NaN NaN NaN NaN NaN NaN
## [4,] -1.311874 NaN NaN NaN NaN NaN NaN NaN NaN
## [5,] -1.331859 NaN NaN NaN NaN NaN NaN NaN NaN
## [6,] -1.352608 NaN NaN NaN NaN NaN NaN NaN NaN
## [,177] [,178] [,179] [,180]
## [1,] NaN NaN NaN NaN
## [2,] NaN NaN NaN NaN
## [3,] NaN NaN NaN NaN
## [4,] NaN NaN NaN NaN
## [5,] NaN NaN NaN NaN
## [6,] NaN NaN NaN NaNavsst = rowMeans(sst,na.rm=FALSE,dims=2)
head(avsst)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] NA NA NA NA NA NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA NA NA NA NA NA
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## [1,] NA NA NA NA NA 3.381521 3.858063 4.336045 4.717195
## [2,] NA NA NA NA 2.640299 3.343127 3.856099 4.347613 4.757762
## [3,] NA NA NA NA 2.605906 3.285104 3.853700 4.337926 4.789467
## [4,] NA NA NA NA NA 3.209619 3.823459 4.288603 4.799861
## [5,] NA NA NA NA NA 3.125613 3.760868 4.249798 4.791644
## [6,] NA NA NA NA NA 3.046399 3.676972 4.190188 4.757562
## [,34] [,35] [,36] [,37] [,38] [,39] [,40]
## [1,] 4.937529 5.121027 5.228368 5.276693 5.383786 5.588450 5.834319
## [2,] 4.992551 5.178441 5.282075 5.332977 5.448545 5.654996 5.893275
## [3,] 5.051107 5.245836 5.352025 5.406577 5.527186 5.732940 5.963246
## [4,] 5.101757 5.308497 5.415242 5.471367 5.595007 5.799395 6.024938
## [5,] 5.134830 5.358847 5.466341 5.523685 5.649658 5.852417 6.075754
## [6,] 5.136606 5.383885 5.499960 5.561589 5.691114 5.894467 6.116872
## [,41] [,42] [,43] [,44] [,45] [,46] [,47]
## [1,] 6.089891 6.392899 6.807126 7.395176 8.198454 9.182922 10.28493
## [2,] 6.142618 6.446611 6.863599 7.448278 8.239780 9.208206 10.29554
## [3,] 6.205951 6.509534 6.930888 7.519525 8.308837 9.270089 10.34655
## [4,] 6.265406 6.569547 6.994052 7.589315 8.382902 9.340843 10.40830
## [5,] 6.317786 6.623745 7.051349 7.656564 8.459444 9.415292 10.47420
## [6,] 6.360271 6.668809 7.099608 7.711499 8.520291 9.473809 10.52468
## [,48] [,49] [,50] [,51] [,52] [,53] [,54]
## [1,] 11.46250 12.66317 13.82411 14.92240 15.94492 16.87619 17.71784
## [2,] 11.46731 12.67373 13.84323 14.94294 15.95125 16.85455 17.67240
## [3,] 11.50894 12.71464 13.88537 14.97584 15.96337 16.84375 17.64730
## [4,] 11.55934 12.75655 13.92118 14.99712 15.96288 16.82730 17.62625
## [5,] 11.61128 12.79043 13.94158 15.00250 15.95163 16.80919 17.61009
## [6,] 11.64938 12.81133 13.94616 14.99220 15.93050 16.78829 17.59514
## [,55] [,56] [,57] [,58] [,59] [,60] [,61]
## [1,] 18.47638 19.17034 19.83273 20.48815 21.13088 21.77575 22.45981
## [2,] 18.42155 19.11552 19.78132 20.44299 21.09587 21.74882 22.43300
## [3,] 18.39475 19.09190 19.75929 20.42428 21.08470 21.74278 22.42178
## [4,] 18.37675 19.07501 19.73970 20.40654 21.07386 21.73468 22.40468
## [5,] 18.36252 19.05641 19.71551 20.38603 21.06221 21.72662 22.38800
## [6,] 18.34791 19.03604 19.69304 20.37081 21.05814 21.72872 22.38244
## [,62] [,63] [,64] [,65] [,66] [,67] [,68]
## [1,] 23.17053 23.84759 24.44503 24.94911 25.36650 25.72611 26.04796
## [2,] 23.14096 23.81549 24.40719 24.90331 25.31559 25.67516 25.99783
## [3,] 23.12099 23.79121 24.38067 24.87531 25.28527 25.64079 25.95887
## [4,] 23.09026 23.75846 24.35744 24.85826 25.26346 25.60874 25.92117
## [5,] 23.05630 23.71929 24.32973 24.83760 25.23590 25.56941 25.87904
## [6,] 23.03004 23.67799 24.29044 24.80547 25.20258 25.52846 25.83635
## [,69] [,70] [,71] [,72] [,73] [,74] [,75]
## [1,] 26.32445 26.55173 26.75948 26.97560 27.17508 27.31499 27.39569
## [2,] 26.27334 26.50220 26.71651 26.93879 27.13534 27.26364 27.33830
## [3,] 26.23155 26.46226 26.68056 26.90190 27.08964 27.20817 27.28132
## [4,] 26.19730 26.43565 26.65797 26.87157 27.04160 27.14863 27.22409
## [5,] 26.16596 26.41685 26.64216 26.84217 26.98924 27.08618 27.16728
## [6,] 26.13248 26.39359 26.61975 26.80794 26.93912 27.03194 27.12083
## [,76] [,77] [,78] [,79] [,80] [,81] [,82]
## [1,] 27.44509 27.50066 27.58782 27.70378 27.83576 27.97417 28.12420
## [2,] 27.39580 27.46339 27.55613 27.67277 27.80587 27.94860 28.10446
## [3,] 27.34564 27.42130 27.52063 27.64253 27.77900 27.92458 28.08102
## [4,] 27.29603 27.37678 27.48196 27.61461 27.76514 27.92178 28.07634
## [5,] 27.24643 27.33084 27.44234 27.58792 27.75747 27.93150 28.08726
## [6,] 27.20686 27.29251 27.40798 27.56726 27.75564 27.94323 28.10084
## [,83] [,84] [,85] [,86] [,87] [,88] [,89]
## [1,] 28.27210 28.39186 28.47415 28.51587 28.52326 28.51095 28.49274
## [2,] 28.25597 28.37095 28.43924 28.46788 28.47116 28.46078 28.44407
## [3,] 28.22859 28.33871 28.40339 28.42512 28.41499 28.38883 28.36379
## [4,] 28.21087 28.30764 28.36446 28.37903 28.35281 28.30614 28.27174
## [5,] 28.20787 28.28541 28.32670 28.33317 28.29583 28.23502 28.19494
## [6,] 28.21430 28.27829 28.30480 28.30233 28.25530 28.18192 28.13556
## [,90] [,91] [,92] [,93] [,94] [,95] [,96]
## [1,] 28.48855 28.52744 28.61702 28.74300 28.88531 29.02458 29.15066
## [2,] 28.44197 28.48622 28.58222 28.71505 28.86327 28.99984 29.11747
## [3,] 28.36995 28.43683 28.55614 28.70197 28.85207 28.98083 29.08622
## [4,] 28.29055 28.38962 28.54268 28.70200 28.84349 28.95996 29.06033
## [5,] 28.22380 28.35119 28.53877 28.71141 28.83930 28.94378 29.04706
## [6,] 28.16775 28.30796 28.51316 28.69661 28.82501 28.92744 29.03195
## [,97] [,98] [,99] [,100] [,101] [,102] [,103]
## [1,] 29.24250 29.27339 29.24505 29.18890 29.11939 29.01377 28.85041
## [2,] 29.20710 29.24266 29.22484 29.18518 29.13362 29.04275 28.88339
## [3,] 29.16990 29.20980 29.20218 29.17470 29.13467 29.05130 28.89371
## [4,] 29.14489 29.19069 29.19297 29.17226 29.13106 29.04408 28.88879
## [5,] 29.13829 29.19097 29.20288 29.18434 29.13087 29.03048 28.87656
## [6,] 29.12959 29.19383 29.21767 29.19998 29.13420 29.02181 28.86612
## [,104] [,105] [,106] [,107] [,108] [,109] [,110]
## [1,] 28.61887 28.31397 27.95036 27.55969 27.17528 26.79081 26.35562
## [2,] 28.64399 28.33023 27.96680 27.58841 27.22257 26.84872 26.40623
## [3,] 28.65287 28.34189 27.98724 27.61537 27.24897 26.86978 26.41777
## [4,] 28.66180 28.37235 28.03292 27.66220 27.28093 26.87706 26.40647
## [5,] 28.67322 28.41692 28.09580 27.72420 27.31910 26.87248 26.37227
## [6,] 28.67588 28.44168 28.13407 27.76081 27.33769 26.86188 26.34389
## [,111] [,112] [,113] [,114] [,115] [,116] [,117]
## [1,] 25.85138 25.31002 24.77544 24.26276 23.74597 23.21039 22.67361
## [2,] 25.88464 25.32946 24.79243 24.28648 23.77252 23.23126 22.68172
## [3,] 25.88461 25.32405 24.78959 24.28694 23.77442 23.22981 22.66252
## [4,] 25.86653 25.31046 24.78290 24.27734 23.75738 23.20474 22.62644
## [5,] 25.82916 25.28718 24.77164 24.26021 23.72834 23.16610 22.58412
## [6,] 25.80253 25.26955 24.75463 24.23170 23.69317 23.13501 22.56103
## [,118] [,119] [,120] [,121] [,122] [,123] [,124]
## [1,] 22.14416 21.59859 21.00727 20.39410 19.82654 19.31977 18.84329
## [2,] 22.13107 21.55917 20.94474 20.32369 19.76499 19.26545 18.77950
## [3,] 22.08592 21.49868 20.88716 20.28231 19.73482 19.22532 18.71058
## [4,] 22.04467 21.46778 20.88108 20.30038 19.75273 19.21278 18.65333
## [5,] 22.01433 21.46285 20.90624 20.34335 19.78324 19.20708 18.60652
## [6,] 22.00482 21.47340 20.93432 20.37059 19.78683 19.18060 18.55815
## [,125] [,126] [,127] [,128] [,129] [,130] [,131]
## [1,] 18.38774 17.95152 17.51797 17.06518 16.60105 16.15039 15.67637
## [2,] 18.30185 17.84551 17.40998 16.98034 16.54860 16.11970 15.66300
## [3,] 18.19802 17.71666 17.27593 16.85351 16.43094 16.01410 15.57471
## [4,] 18.10097 17.59691 17.14795 16.72215 16.29981 15.89102 15.46969
## [5,] 18.02028 17.49466 17.02976 16.59160 16.16698 15.76438 15.35857
## [6,] 17.95510 17.41028 16.92409 16.47031 16.04095 15.64127 15.24576
## [,132] [,133] [,134] [,135] [,136] [,137] [,138]
## [1,] 15.07315 14.25427 13.27339 12.30132 11.46592 10.75858 10.08537
## [2,] 15.08177 14.29197 13.34771 12.40458 11.58651 10.89705 10.24591
## [3,] 15.02137 14.27896 13.39439 12.50756 11.73678 11.08742 10.47371
## [4,] 14.94928 14.26540 13.45987 12.64642 11.93071 11.32665 10.75543
## [5,] 14.87204 14.25147 13.53148 12.79128 12.12069 11.54936 11.01216
## [6,] 14.78554 14.21496 13.56059 12.87762 12.23904 11.68490 11.17282
## [,139] [,140] [,141] [,142] [,143] [,144] [,145]
## [1,] 9.416843 8.891008 8.423951 8.030039 7.617577 7.159713 6.653107
## [2,] 9.596250 9.073308 8.596709 8.180797 7.740390 7.254704 6.721100
## [3,] 9.851744 9.327616 8.824090 8.361171 7.873391 7.350591 6.787400
## [4,] 10.165143 9.634758 9.085650 8.554732 8.009507 7.441781 6.839059
## [5,] 10.457214 9.929723 9.330962 8.727171 8.125147 7.512736 6.867696
## [6,] 10.651357 10.134795 9.509652 8.851994 8.201692 7.549889 6.870418
## [,146] [,147] [,148] [,149] [,150] [,151] [,152]
## [1,] 6.096617 5.477840 4.792800 4.051160 3.288124 2.563124 1.893136
## [2,] 6.135208 5.488662 4.782238 4.025948 3.254650 2.524728 1.874019
## [3,] 6.176583 5.504535 4.772963 3.998688 3.217939 2.484215 1.842973
## [4,] 6.197486 5.503200 4.754645 3.970628 3.186019 2.449848 1.811585
## [5,] 6.189107 5.474352 4.721331 3.938311 3.154519 2.412405 1.770459
## [6,] 6.158705 5.422975 4.665592 3.886683 3.107264 2.364477 1.726514
## [,153] [,154] [,155] [,156] [,157] [,158] [,159] [,160]
## [1,] 1.308794 0.7784685 0.2631246 NA NA NA NA NA
## [2,] 1.306799 0.7908744 0.3058383 -0.2135148 NA NA NA NA
## [3,] 1.285715 0.7773835 0.3033367 -0.2135408 NA NA NA NA
## [4,] 1.254575 0.7492272 0.2851924 NA NA NA NA NA
## [5,] 1.210653 0.7045036 0.2491880 NA NA NA NA NA
## [6,] 1.158537 0.6669740 0.2311133 -0.2632437 NA NA NA NA
## [,161] [,162] [,163] [,164] [,165] [,166] [,167] [,168] [,169]
## [1,] NA NA NA NA NA NA NA -1.232284 NA
## [2,] NA NA NA NA NA NA NA -1.264675 NA
## [3,] NA NA NA NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA NA NA NA
## [,170] [,171] [,172] [,173] [,174] [,175] [,176] [,177] [,178] [,179]
## [1,] NA NA NA NA NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA NA NA NA NA
## [,180]
## [1,] NA
## [2,] NA
## [3,] NA
## [4,] NA
## [5,] NA
## [6,] NAFrom the output it is clear that na.rm = TRUE removes all the NA’s from the dataframe.
Two lines of code to get nicer colours for the eventual plots
colors <- rev(brewer.pal(10, "RdYlBu"))
pal <- colorRampPalette(colors)for this we can use ‘levelplot’ in the lattice package simple e.g.
levelplot(avsst,col.regions = pal(100));
instantly we have a colorbar and the plot looks nice however, adding lat and lon is a little more involved as levelplot likes to work with data frames
# expand.grid will get lats and lons added to a data frame expanded on each others size
grid <- expand.grid(x=lon, y=lat)
# our avsst that we want to plot can be added as a vector, basically expanding what we had
grid$avsst <- as.vector(avsst)
levelplot(avsst~x*y,grid,col.regions = pal(100),
xlab='Longitude',ylab='Latitude',main='Average SST'
) +
layer(sp.lines(world.coast))
From the output it is clear that SST is much warmer along the equator in comparison to poleward regions. This may be explained by circulation in the ocean.
We’re not interested in short term variability so let’s make annual averages
yrs <- annmean$Group.1
nyr <- length(yrs)
asst <- array(NA,c(dim(lon),dim(lat),nyr)) # asst will contain the annual sst, declare it here to have the correct size, this is to speed up computation
# use a for loop to load the annual data into asst
# note: for loops are a bit old school and I'm sure there is a neat vectorised way to do thisâif you know it, let me know
for (k in 1:nyr) {
asst[,,k] <- rowMeans(sst[,,year==yrs[k]],na.rm=FALSE,dims=2) # annual averages from monthly data
}note that all we’ve done is create annual averages so the average looks exactly the same add to the same data frame as earlier for plotting
grid$an_avsst <- as.vector(rowMeans(asst,na.rm=FALSE,dims=2))
#plot
levelplot(an_avsst~x*y, data=grid,col.regions = pal(100),xlab='Longitude',ylab='Latitude',main='Annually Averaged SST') +layer(sp.lines(world.coast)) 
From the annually averaged SST it is clear that it is the warmest along the equator, while the poles are the coldest.
remove the global mean from each year lets do a traditional loop
gmean <- colMeans(asst, na.rm = TRUE, dims=2)
for (k in 1:nyr){
asst[,,k]<-asst[,,k]-matrix(gmean[k],length(lon),length(lat))
}in this example, we’re going to take a timeseries of sst at a point lon0, lat0 and look at the patterns of correlation of sst with this tseries will just be the asst value at a given point lon0, lat0 near ireland lon0 <- -10.5, lat0 <- 51.5 subpolar gyre lon0 <- -35.5, lat0 <- 59.5 try a few of your own note: keep the .5 at the end to find a matching lat/lon
lon0 <- -10.5 #
lat0 <- 51.5 #
sst_ts<-asst[which(lon==lon0),which(lat==lat0),]look at this timeseries:
plot(yrs,sst_ts,type='l',xlab='Year',ylab='SST Anomaly',main=paste0('SSTA at Long=', lon0, ',Lat=', lat0))
lets calculate the correlation in a loop also can you remove one of these loops? does this speed up computation time? hints: help(cor)
cmatrix <- matrix(NA,dim(lon),dim(lat))
for (i in 1:dim(lon)) {
for (j in 1:dim(lat)) {
cmatrix[i,j] <- cor(asst[i,j,], sst_ts)
}
}add to the same data frame as earlier for plotting
grid$corr <- as.vector(cmatrix)plot
levelplot(corr~x*y, data=grid , xlim=c(-120,10),ylim=c(0,80), # at=c(-1:1),
col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SSTA with Long=', lon0, ',Lat=', lat0)) +
layer(sp.lines(world.coast)) 
task 2
repeat the correlation plot using a point north of the gulf stream
lon0 <- -60.5 #
lat0 <- 39.5 #
sst_ts<-asst[which(lon==lon0),which(lat==lat0),]
cmatrix <- matrix(NA,dim(lon),dim(lat))
for (i in 1:dim(lon)) {
for (j in 1:dim(lat)) {
cmatrix[i,j] <- cor(asst[i,j,], sst_ts)
}
}
grid$corr <- as.vector(cmatrix)
levelplot(corr~x*y, data=grid , xlim=c(-120,10),ylim=c(0,80), # at=c(-1:1),
col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('Correlation of SSTA with Long=', lon0, ',Lat=', lat0)) +
layer(sp.lines(world.coast)) 
Using a grid point north of the Gulf stream changes the map completely. It appears much cooler than using a grid point below the gulf stream.
task 3: rather than using “cor” repeat using “lm”. Calculate the slope at each point and plot the pattern
lon0 <- -60.5 #
lat0 <- 39.5 #
sst_ts<-asst[which(lon==lon0),which(lat==lat0),]
r.matrix <- matrix(NA, length(lon), length(lat))
for (i in 1:dim(lon)) {
for (j in 1:dim(lat)){
if(length(asst[i,j,][!is.na (asst[i,j,])])>2){
r.lm <- lm(asst[i,j,]~sst_ts)
r.matrix[i,j] <- r.lm$coefficients[2]
}
}
}
grid$lm <- as.vector(r.matrix)
levelplot(lm~x*y, data=grid , xlim=c(-120,10),ylim=c(0,80), at=c(-20:20)/10,
col.regions = pal(100),xlab='Longitude',ylab='Latitude',main=paste0('LM of SSTA with Long=', lon0, ',Lat=', lat0)) +
layer(sp.lines(world.coast)) 