-
Notifications
You must be signed in to change notification settings - Fork 0
/
loadncdf.r
72 lines (49 loc) · 1.95 KB
/
loadncdf.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
# convert netcdf to raster http://geog.uoregon.edu/GeogR/topics/netcdf-to-raster.html
require(ncdf4)
library(raster)
library(maptools)
library(maps)
library(rasterVis)
library(dplyr)
workdir <-'//Users//mathewbrown//projects//Rstuff//carbonFlux'
setwd(workdir)
cname <- 'CT2016.flux1x1.200101.nc'
ncin <- nc_open(cname)
print(ncin)
#sum flux components for total surface co2 exchange
#tot = ncatt_get(ncin,'bio_flux_opt')
# use raster to read ncdf
bio <- raster(cname,varname='bio_flux_opt')
ocn <- raster(cname,varname='ocn_flux_opt')
fossil <- raster(cname,varname='fossil_flux_imp')
fire <- raster(cname,varname='fire_flux_imp')
monthly=12*3600*24*365
#Put data into data frame
mapPointsBio <- as.data.frame(rasterToPoints(bio))
mapPointsOcn <- as.data.frame(rasterToPoints(ocn))
mapPointsFossil <- as.data.frame(rasterToPoints(fossil))
mapPointsFire <- as.data.frame(rasterToPoints(fire))
#merge data frames
all=full_join(mapPointsBio,mapPointsOcn)
all=full_join(all,mapPointsFossil)
all=full_join(all,mapPointsFire)
#jet get biosphere flux
mapPointsBio$bio_flux_opt = mapPointsBio$bio_flux_opt*monthly
write.csv(mapPointsBio, file=paste0(substr(cname,1,nchar(cname)-3),".csv"), row.names=F)
df_sum= all %>%
rowwise() %>%
mutate(tot=sum(bio_flux_opt,ocn_flux_opt))
df_sum$bio_flux_opt = df_sum$bio_flux_opt*monthly
df_sum$tot = df_sum$tot*monthly
write.csv(df_sum, file=paste0(substr(cname,1,nchar(cname)-3),".csv"), row.names=F)
datain <- raster(cname,varname='bio_flux_opt')
datain <- setMinMax(datain)
print(datain)
saveRDS(datain,file=paste0("200101",'.rds'))
# plot an example
world.outlines <- map('world',plot=FALSE)
world.outlines.sp <- map2SpatialLines(world.outlines, proj4string = CRS("+proj=longlat"))
mapTheme <- rasterTheme(region = rev(brewer.pal(10, "RdBu")))
plt <- levelplot(datain, margin = F, cuts=30, pretty=TRUE, par.settings = mapTheme,
main="Carbon Flux")
plt + layer(sp.lines(world.outlines.sp, col = "black", lwd = 0.5))