Suite

Vous voulez des valeurs de régression linéaire de cellule pour un netCDF ou un raster multi-canaux


J'ai des données netCDF de type raster climatique.

Je veux la pente et d'autres statistiques de régression de la température pour chaque cellule au fil du temps.

J'ai lu que GEOV IDL peut faire une régression temporelle pour obtenir ce que je veux mais ce n'est pas open source. Est-ce que quelqu'un connaît un bon logiciel open source qui peut accomplir ce que je voudrais ? R peut-être ?


Voici comment vous pouvez obtenir la pente, en utilisant R et le package raster. Pour (aussi) obtenir l'interception voir help(calc)

library(raster) # votre fichier # b <- brick("file.nc") # exemple de données : b <- brick(system.file("external/rlogo.grd", package="raster")) # here time est de 1 à n, mais vous pouvez le définir autrement time <- 1:nlayers(b) # écrivez une fonction qui réexécute la ou les valeurs d'intérêt fun <- function(x) { lm(x ~ time)$coefficients[ 2] } x <- calc(b, fun) plot(x)

Si vous avez des valeurs NA, vous avez besoin d'une fonction plus complexe

fun2 <- function(x) { d <- na.omit(cbind(x, time)) if (nrow(d) > 2) { lm(x ~ time, data=data.frame(d))$coefficients[ 2] } else { NA } } b[1:10] <- NA x2 <- calc(b, fun2)

pour obtenir r^2

amusant <- function(x) { summary(lm(x ~ time))$r.squared }

Addenda

J'ai remarqué que cette fonction est un peu lente car, pour chaque cellule raster, elle correspond à un modèle et renvoie beaucoup d'informations qui ne sont pas utilisées. Pour chaque modèle, certains calculs sont répétés (car la variable indépendante est fixe). Si vous voulez une sortie simple comme la pente ou l'interception, vous pouvez facilement raccourcir les choses en les calculant directement via l'algèbre linéaire et en pré-calculant certains des résultats intermédiaires (constants).

Pour le cas sans NA uniquement :

library(raster) b <- brick(system.file("external/rlogo.grd", package="raster")) time <- 1:nlayers(b) LMfun <- function(x) { lm(x ~ time )$coefficients[2] } system.time( xlm <- calc(b, LMfun) ) # user system elapsed # 7.95 0.00 7.96 # ajouter 1 pour un modèle avec une interception X <- cbind(1, time) # pre- calcul de la partie constante des moindres carrés invXtX <- solve(t(X) %*% X) %*% t(X) # modèle de régression très réduit. [2] est d'obtenir la pente LAfun <- function(y) (invXtX %*% y)[2] system.time( xla <- calc(b, LAfun) ) # user system elapsed # 0.06 0.00 0.06

Cette approche est donc environ 130 fois plus rapide !


Si vous êtes familier avec Python, vous pouvez utiliser la bibliothèque netCDF4-python qui peut lire et écrire les données netCDF 3 et 4 surnumpytableaux. Par exemple:

from netCDF4 import Dataset root_group = Dataset("path_to_dataset", format="NETCDF4") print root_group $ netCDF4 style dump data = root_group.variables["some_variable"][:]

Python possède un grand nombre d'autres bibliothèques que vous pouvez utiliser pour la régression. Personnellement, mon premier arrêt serait scikit-learn, mais vous devriez également essayer les modèles de statistiques (surtout si vous connaissez R - les modèles sont définis de manière très similaire).

Je ne suis pas aussi familier avec R, mais il existe certainement un package ncdf4 de CRAN et vous pouvez utiliser le lm méthode pour ajuster facilement des modèles linéaires par exemple. De nombreux tutoriels sont également disponibles.


GDAL prend en charge les fichiers netCDF :

gdalinfo --formats|grep -i cdf GMT (rw): GMT NetCDF Grid Format netCDF (rw+): Network Common Data Format

Ainsi, vous pouvez ouvrir ce genre de fichiers directement dans QGIS. Pour cette page :

http://www.unidata.ucar.edu/software/netcdf/examples/files.html

J'ai téléchargé cet exemple de fichier climatique ECMWF_ERA-40_subset.nc. Il a 17 couches multibandes (chaque couche avec 62 bandes : 1054 au total) et elles ont été totalement chargées dans QGIS. J'ai rendu la couche p2t ("température de 2 mètres") en pseudo-couleur (5 classes) et le résultat était :

On peut observer que l'image est centrée dans l'océan Pacifique avec un rectangle délimité de : 0 90 357,5 -90 (Voir le fichier de métadonnées). Il s'agit d'une donnée importante pour attribuer une projection en utilisant gdalwarp dans QGIS (Raster -> Projections -> Warp).

Importé dans GRASS, vous pouvez calculer la régression linéaire à partir de deux cartes raster (y = a + b*x) avec r.regression.line.

D'autre part, avec l'aide de Value Tool Plugin dans QGIS, vous pouvez effectuer une analyse précédente exploratoire de tous les jeux de données et bandes :