1. install.packages("raster")
2. install.packages("ncdf4")
3. install.packages("fields")
4. install.packages("sp")
5. install.packages("rasterVis) # raster visualisation

1. install.packages("plotly")
2. install.packages("ggplot2")
3. install.packages("maptools")
4. install.packages("rgeos")
5. install.packages("geoR")
6. install.packages("maps")

library(ncdf4) library(fields) library(sp) # classes for spatial data library(raster)

library(rasterVis) # raster visualisation library(plotly) library(ggplot2) library(maptools) library(rgeos) library(geoR) library(maps)

data(wrld_simpl)

titula="LGM annual precipitation, UKMO simulation" month=1

1. 1. urli="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/gcmoutput/pmip/ukmo/21cal/ukmo_21cal_pr_ca_map.nc"

urli="https://www.ncei.noaa.gov/pub/data/paleo/gcmoutput/pmip/ukmo/21cal/ukmo_21cal_pr_ca_map.nc"'

fname <- "/home/mikkonen/r0/ukmo_21cal_pr_ca_map.nc"

download.file(urli, fname, method='auto', quiet = FALSE, mode = "w",

             cacheOK = TRUE,
             extra = getOption("download.file.extra"))

ncfilee <- nc_open(fname) data <- ncvar_get(ncfilee)

longitude <- ncvar_get(ncfilee,"longitude") latitude <- ncvar_get(ncfilee,"latitude")

latitude <-rev(latitude)

pr<-ncvar_get(ncfilee,"pr")

pr<-pr[,c(ncol(pr):1)]

pr365=pr*365

lon=longitude

k=which(lon==180) px=c(k:length(lon),1:(k-1)) lon2=lon[px] pr3652o=pr365[px,,] lon2[1:(which(lon2==0.0))-1]=lon2[1:(which(lon2==0.0)-1)]-360

xloppu=nrow(pr) xalku=1 xkeski1=xloppu/2 xkeski2=xkeski1+1

d=pr365[xalku:xkeski1,]

e=pr365[xkeski2:xloppu,]

pr3652=rbind(e,d)

paletti=rev(topo.colors(12))

tasot<-c(0,100,200,400,600,800,1200,1600,2400,3200,6400)

omapaletti<-c("#5D00FF", "#002EFF","#00B9FF","#00FFB9" ,"#00FF2E","#5DFF00","#E8FF00", "#FF8B00","red", "#FF008B","#E800FF")

png(file="/home/mikkonen/r0/myplot.png", width=1200, height=800)

filled.contour(lon2, latitude, pr3652,

key.title = title(main = "mm"),

key.axes = axis(4, tasot ), plot.title=title(main=titula,xlab="longitude", ylab="latitude"), col=paletti, levels=tasot, plot.axes = { contour(lon2, latitude, pr3652, levels=tasot, nlevels =12 ,

drawlabels = TRUE, axes = FALSE, frame.plot = FALSE, add = TRUE);

axis(1, seq(-180, 180, by = 30)); axis(2, seq(-90, 90, by = 30));

map('world', add=TRUE); } )

dev.off()

1. 1. process ukmo pmpip 2003 output precip file

library(raster)

iname0<-"./data1/ukmo_0cal_pr_ca_map.nc"

1. iname0<-"./data1/ukmo_21cal_pr_ca_map.nc"

precip0<-raster(iname0)

print(precip0)

raster::plot(precip0)

precip1<-rotate(precip0) precip1=precip1*365

crs(precip1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

ext1<-c(-10,10,40,50)

precip2<-crop(precip1, ext1)

raster::plot(precip1)

crs(precip2) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

outname1 = "precip1.nc" outname2 = "precip2.nc"

raster::writeRaster(precip1, outname1, format="CDF", varname="pr", overwrite= TRUE)

raster::writeRaster(precip2, outname2, format="CDF", varname="pr", overwrite= TRUE)