# na poczatek tworzymy polaczenie ze zrodlowym plikiem NetCDF'a
# ktory zawiera koordynaty domeny WRF'a
setwd("/home/bartosz/Pulpit/natalia/")
library(ncdf4) # biblioteka do NetCDF'ow
nc <- nc_open("us22.nc") # na poczatek wyciagamy pierwszy, jakis obliczony juz plik
zmienna <- ncvar_get(nc=nc,varid="U10") #wczytujemy jakas wyliczona zmienna do obiektu
nc_close(nc) # zamykamy polacenie do NetCDF'a
dim(zmienna) # sprawdzmy jeszcze wymiary naszego pliku...
## [1] 276 276 24
image(zmienna[,,1]) # wyglada ok...
# problem z tym plikiem polega na tym, ze nie ma w nim koordynat!!!
##############
##############
# sprobujmy zatem wczytac dowolny plik WRFOUT w ktorym te dane sa...
nc <- nc_open("wrfout_d03_2009-01-31_00:00:00") # na poczatek wyciagamy pierwszy, jakis obliczony juz plik
# a nastepnie wyciagniemy z tego pliku zmienne dla dlugosci i szerokosci geograficznej
lon <- ncvar_get(nc=nc,varid="XLONG")
lat <- ncvar_get(nc=nc,varid="XLAT")
# i tutaj sie pojawia pierwsza niespodzianka...
dim(lon)
## [1] 276 276 24
# normalnie dla rastra wystarcza podpisy wartosci w formie jednowymiarowego wektora (jako podpis osi)
# ... tutaj dostajemy macierz wartosci, poniewaz siatka WRFa jest nieregularna!!!
identical(lon[,,1], lon[,,2]) # wszystkie wymiary maja takie same wartosci, zatem wystarcza nam
## [1] TRUE
# 2 wymiary zamiast 3:
lon <- lon[,,1]
lat <- lat[,,1]
library(fields)
## Loading required package: spam
## Loading required package: grid
## Spam version 1.4-0 (2016-08-29) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
## Loading required package: maps
image.plot(lon) # w ten sposob zmieniaja sie wartosci poludnikow
image.plot(lat) # a tak rownoleznikow
## w najprostszej postaci mozemy przekonwertowac nasze dane do ramki danych w formacie X-Y-Z
xy1 <- data.frame(lon=as.numeric(lon), lat=as.numeric(lat), zmienna=as.numeric(zmienna))
# musimy zmienic klase obiektu na bardziej 'GIS'owa:
# potrzebne nam beda biblioteki:
library(sp)
library(rgdal)
## rgdal: version: 1.2-8, (SVN revision 663)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
## Path to GDAL shared files: /usr/share/gdal/2.2
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.2-5
coordinates(xy1) <- ~lon+lat
proj4string(xy1) <- CRS("+init=epsg:4326") # nadajemy projekcje WGS-84 (czyli klasyczne lon-lat)
head(xy1)
## coordinates zmienna
## 1 (3.309891, 74.40241) -1.413190
## 2 (3.405884, 74.40943) -1.335547
## 3 (3.501968, 74.41634) -1.263924
## 4 (3.598129, 74.42326) -1.240898
## 5 (3.694382, 74.43009) -1.266324
## 6 (3.790718, 74.43693) -1.289069
spplot(xy1["zmienna"], scales=list(draw=T), main="tytul") # to moze troche czasu zajac
# bo teraz rysujemy kazdy z 76tys. punktow oddzielnie
# OK., Ja nie przepadam za projekcja WGS 84 w obszarach polarnych i wole np epgs:3575
# Transformacje wspolrzednych mozna wykonac za pomoca:
xy1_epsg3575 <- spTransform(xy1, CRS("+init=epsg:3575"))
spplot(xy1_epsg3575) # i jeszcze raz sobie rysujemy w nowym ukladzie wspolrzednych
# wszystko co do tej pory rysowalismy bazowalo na formacie wektorowym, podczas
# kiedy najprawdopodobniej chcialabys uzyskac regularna siatke rastra
library(rasterVis)
## Loading required package: raster
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
library(raster)
naszraster <- raster(xmn=extent(xy1_epsg3575)[1], xmx=extent(xy1_epsg3575)[2] ,
ymn=extent(xy1_epsg3575)[3] , ymx=extent(xy1_epsg3575)[4],
ncols=276, nrow=276,crs=CRS("+init=epsg:3575")) # tworzymy raster o naszych wymiarach i projekcji
naszraster <- rasterize(xy1_epsg3575, naszraster, field="zmienna", fun=mean, method='bilinear')
image(naszraster) # voila...
# lub:
plot(naszraster)
# gdybyś jednak chciała pozostać przy WGS-84 i uzyskać raster w WGS-84
naszraster_wgs84 <- raster(xmn=extent(xy1)[1], xmx=extent(xy1)[2] ,
ymn=extent(xy1)[3] , ymx=extent(xy1)[4],
ncols=276, nrow=276,crs=CRS("+init=epsg:4326")) # tworzymy raster o naszych wymiarach i projekcji
naszraster_wgs84 <- rasterize(xy1, naszraster_wgs84, field="zmienna", fun=mean, method='bilinear')
image(naszraster_wgs84) # voila...
library(mapdata)
map('world',add=T) # dodajemy kontury
# a jesli mamy wlasne pliki SHP z konturami to mozna to rozwiazac np. tak:
swiat <- shapefile("mapa_svalbard.shp")
proj4string(swiat) <- CRS("+init=epsg:4326") # zakladamy ze standardowo nasz zbior danych ma projekcje WGS-84
## Warning in `proj4string<-`(`*tmp*`, value = <S4 object of class structure("CRS", package = "sp")>): A new CRS was assigned to an object with an existing CRS:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## without reprojecting.
## For reprojection, use function spTransform
swiat <- spTransform(swiat, CRS("+init=epsg:3575")) # zmieniamy uklad z WGS-84 na epsg:3575
plot(naszraster, col=rainbow(64))
plot(swiat, add=T, lwd=3)
