search

mdsumner / python-aos-lesson

Python (and git) lesson for atmosphere and ocean scientists
https://carpentrieslab.github.io/python-aos-lesson/
Other
0 stars 0 forks source link

Idioms from Python #1

Open mdsumner opened 5 years ago

mdsumner commented 5 years ago

Some examples from python we could use in R

https://carpentrieslab.github.io/python-aos-lesson/

xarray has this for a climatology over time

dset['pr'].mean('time')

In R with a 3D array (x, y, time) that would be

apply(pr, 1:2, fun = mean)

and in raster it's

calc(pr, mean)
## or
stackApply(pr, mean)

I'm not sure what raster needs to keep it out of memory, but at any rate we can't set margin to anything but 1:2 in raster.

group by

Here "time.season" appears to be an xarray thing on an intepreted time axis, it's the triplets of months? (That seems common, but I don't see these entities in the ACCESS files).

dset['pr'].groupby('time.season').mean('time')

In R, for the particular file 

time <- as.Date(getZ(dset))
time_season <- function(x) {
  setNames(rep(c("DJF", "MAM", "JJA", "SON"), each = 3), month.abb[c(12, 1:11)])[format(x, "%b")]
}

## we have to make sure the codes align, they won't here - perhaps use ordered
 stackApply(dset, as.integer(factor(time_season(as.Date(getZ(dset))))), fun = mean)
mdsumner commented 5 years ago

There's a lot of this in the package remote

mdsumner commented 5 years ago

Code to reproduce the first plot, (close)

access_pr_file <- 'data/pr_Amon_ACCESS1-3_historical_r1i1p1_200101-200512.nc'

library(raster)
dset <- brick(access_pr_file, varname = "pr")
clim <- stackApply(dset, indices = rep(1, nlayers(dset)), fun = mean) * 24 * 3600

## the default plotting in python seems to clamp to this range, based on the
## input breaks
clim[clim > 12.5] <- 12.5
plot(raster(list(x = lon, y = lat, z = x)), col = viridis::viridis(11), zlim = c(0, 12.5))
maps::map("world2", add = TRUE)