mariusz_kupczyk_gisR.R

#######################################################################################################################
library(rworldmap)

## Loading required package: sp

## ### Welcome to rworldmap ###

## For a short introduction type :   vignette('rworldmap')

library(ncdf4)
##########################################################PACZKI#######################################################

newmap <- getMap(resolution = "low")

#####################################################

################ZDEFINIOWANIE PLIKU#################
ncin <- nc_open('~/SBCAPE.nc')
###################################################



######################################DLUGOSC I SZEROOSC###############################
lon <- ncvar_get(ncin, "lon")
lat <- ncvar_get(ncin, "lat", verbose = F)
dname="SBCAPE"
dunits <- ncatt_get(ncin, dname, "units")
t <- ncvar_get(ncin, "time")
tunits <- ncatt_get(ncin, "time", "units")
########################################################################################


tmp.array <- ncvar_get(ncin, "SBCAPE")
czas <- as.POSIXct(t, origin="1970-01-01", tz='GMT')


#maxm <- dim(czas)
m <- 10
#for (m in 1550:1600){
  
  tmp.slice <- tmp.array[, , m]
  grid <- expand.grid(lon = lon, lat = lat)
  
  
  #png(filename=paste("pliki/d",m,".png", sep=""), width = 900, height = 650) 
  
  colfunc <- colorRampPalette(c("white","cyan","green","yellow","orange", "red", "#600000"))
  
  filled.contour(lon, sort(lat), tmp.slice, color.palette=colfunc, zlim=c(0,5500), ylab="latitude", xlab="longitude",
                 plot.axes = {plot(newmap,
 key.title = title(main = paste(format(czas[m],"%Y-%m-%d %H:%M"),paste("[",dunits$value,"]", sep=""))),
 add=T, border="grey")+
      contour(nlevels=3,labcex=0.85,sort(lon), sort(lat), tmp.slice, add=T, lwd=0.5, col = par("fg"))})

#  dev.off()
#  graphics.off()
  print(format(czas[m],"%Y-%m-%d %H:%M"))

## [1] "1979-01-03 06:00"

#}

################################################################################

lonlat <- expand.grid(lon, sort(lat))
tmp.vec <- as.vector(tmp.slice)
lonlat$z <- tmp.vec

coordinates(lonlat) <- ~Var1+Var2

library(rasterVis)

## Loading required package: raster

## Loading required package: lattice

## Loading required package: latticeExtra

## Loading required package: RColorBrewer

## proba ze stworzeniem rastra:
X=lon
Y=lat
r <- raster(ncol=length(X)-1, nrow=length(Y)-1, xmn=X[1], 
            xmx=X[length(X)], ymn=Y[length(Y)], ymx=Y[1]) ## stworzenie obiektu rastrowego
r

## class       : RasterLayer 
## dimensions  : 58, 107, 6206  (nrow, ncol, ncell)
## resolution  : 0.75, 0.75  (x, y)
## extent      : -27, 53.25, 31.5, 75  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

#res(r) <- 0.75 # nadanie rozdzielczosci dla siatki

library(maptools)

## Checking rgeos availability: FALSE
##      Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
##      which has a restricted licence. It is disabled by default;
##      to enable gpclib, type gpclibPermit()

europa <- readShapeSpatial("/home/bartosz/Dropbox/Dokumenty/projekty/isok/shp/europa.shp")

## Warning: use rgdal::readOGR or sf::st_read

## Warning: use rgdal::readOGR or sf::st_read

wgs84 = CRS('+init=epsg:4326')
#epsg2180 <- CRS('+init=epsg:2180')
#spplot(pl)

proj4string(lonlat) <- wgs84 # na wszelki wypadek gdyby nie bylo jasne...
proj4string(europa) <- wgs84 # j.w.
#wynik33 <- spTransform(lonlat, CRS("+init=epsg:2180")) # mozna ewentualnie transformowac do innego ukladu wspolrzednych:

warstwa <- rasterize(lonlat, r)
warstwa <- warstwa[[2]]
proj4string(warstwa) <- wgs84

# rozpoczecie rysowania
p <- spplot(warstwa, at=c(0,100,200,400,600), col.regions=c("white","yellow","lightgreen","green"),
            scales=list(draw=TRUE), xlab="długość geograficzna", main="PRZYKŁADOWY TYTUŁ MAPY",
            ylab="SZEROKOŚĆ GEOGRAFICZNA")
p <- p + layer(sp.polygons(europa, lwd=0.8, col='gray',fill=F))
p <- update(p+contourplot(warstwa, contour=T, 
                          at=c(200,400,600), labels=T,lwd=0.8))
print(p)