utilisation de "colorbar.plot" commune pour les c

[Quynh-Mai Chu]

Bonjour,

Je bloque à insérer une légende commune pour 4 cartes climatiques. Voici l'image sur le lien suivant:

https://dl.dropboxusercontent.com/u/529 ... d_help.png

Voici un morceau de code produisant l'image:

Code : Tout sélectionner


library(grDevices);
library(clim.pact);
library(fields);

couleur <- c(rgb(255/255,255/255,178/255), rgb(254/255,204/255,92/255), rgb(253/255,141/255,60/255), rgb(227/255,26/255,28/255));
color_plate <- colorRampPalette(couleur, space="rgb");
brks <- seq(min(tas), max(tas), length.out=100);

options(device=X11());
par(mfrow=c(2,2), oma=c(0,0,0,3));

for(i in seq(1,4)){
 
  image(lon_transf, lat, tas_transf[,,i+10*(i-1)], col=color_plate(length(brks)-1), breaks=brks);
  addland(col="black", lwd=1.5);
 
}

colorbar.plot(x=100, y=10, strip=brks, horizontal=T, col=color_plate(length(brks)-1), legend.only=T, lab.breaks=brks[seq(1,100,10)]);



La légende est située en plein milieu de la 4ème carte. Je voudrais faire plusieurs choses:

* Essayer de placer la légende à l'extérieur de la 4ème carte pour en faire une légende commune. Elle disparait lorsqu'elle se situe hors de la carte.

* Agrandir la taille de la légende qui prend toute la bordure droite de l'image

* Ajouter les valeurs à la barre de légende (je pensais que c'était l'argument "legend.breaks" mais rien ne s'affiche)

* Mettre un symbole du genre "°C" en haut à droite de la légende.

Toute idée est la bienvenue :)

Merci beaucoup :)

[Renaud Lancelot]

Il nous faudrait des données pour reproduire l'exemple. S'il s'agit de matrices de données, j'utiliserais les fonctions du package raster pour gérer les données (stack etc), puis plot ou spplot pour les représenter.

[Quynh-Mai Chu]

Bonjour Renaud,

Voici le fichier de données:

https://dl.dropboxusercontent.com/u/529 ... -200911.nc

Elle est sous format netcdf. Je vais mettre le code en entier:

Code : Tout sélectionner


# setwd("~/Dropbox/Master1/ComplModClim/Travail/Data");
setwd("G:/Drop/Dropbox/Master1/ComplModClim/Travail/Data")

library(ncdf); # pour les donnees en format netcdf

data <- open.ncdf("tas_Amon_EC-EARTH_historical_r1i1p1_194001-200911.nc");

t <- get.var.ncdf(data, "time"); # Extraction des variables
lon <- get.var.ncdf(data, "lon");
lat <- get.var.ncdf(data, "lat");
tas <- get.var.ncdf(data, "tas");

nlon <- length(lon);
n <- 1:nlon;
nlat <- length(lat);
Nt <- length(t);

# Transformation des coordonnees longitudes en -180 a 180

lon1 <- lon[n > nlon/2]-360;
lon2 <- lon[n <= nlon/2];
lon_transf <- c(lon1, lon2);

tas_transf <- array(dim=c(nlon, nlat, Nt));
tas_transf[n <= nlon/2,,] <- tas[n > nlon/2,,];
tas_transf[n > nlon/2,,] <- tas[n <= nlon/2,,];

library(grDevices); # pour la fonction colorRampPalette
library(clim.pact); # pour la fonction addland
library(fields); # pour la fonction image


couleur <- c(rgb(255/255,255/255,178/255), rgb(254/255,204/255,92/255), rgb(253/255,141/255,60/255), rgb(227/255,26/255,28/255));
color_plate <- colorRampPalette(couleur, space="rgb");
brks <- seq(min(tas), max(tas), length.out=100);

options(device=X11());
par(mfrow=c(2,2), oma=c(0,0,0,3));

for(i in seq(1,4)){
 
  image(lon_transf, lat, tas_transf[,,i+10*(i-1)], col=color_plate(length(brks)-1), breaks=brks);
  addland(col="black", lwd=1.5);
 
}

# C'est cette commande qui pose problème grr :p

colorbar.plot(x=50, y=10, strip=brks, horizontal=F, col=color_plate(length(brks)-1), legend.only=T, lab.breaks=as.character(brks[seq(1,100,10)]), adj.x=0.5, adj.y=0.5);



Si j'ai vraiment le choix, je préfère pas utiliser le package raster. Il m'a l'air tellement compliqué pour ce que je fais. Je ne fais que des cartes climatiques, pas de la télédétection. A moins que tu connais un site facile qui explique bien l'utilisation du paquet et qui concerne mon cas, je suis preneuse :-).

[Quynh-Mai Chu]

Finalement, je me suis mise à l'utilisation du paquet "raster" (malgré les débuts difficiles, Le résultat est plutôt satisfaisant :-)). J'obtiens bien ce que je veux concernant la légende. Cependant, un nouveau problème surgit: je n'arrive pas à ajouter les continents.

Voici l'image:

https://dl.dropboxusercontent.com/u/529 ... raster.png

Voici le code:

Code : Tout sélectionner



setwd("G:/Drop/Dropbox/Master1/ComplModClim/Travail/Data");

library(ncdf);

data <- open.ncdf("tas_Amon_EC-EARTH_historical_r1i1p1_194001-200911.nc");
lon <- get.var.ncdf(data, "lon");
lat <- get.var.ncdf(data, "lat");
tas <- get.var.ncdf(data, "tas");
t <- get.var.ncdf(data, "time");

nlon <- length(lon);
nlat <- length(lat);
n <- 1:nlon;

Nt <- length(t);

lon1 <- lon[n > nlon/2]-360;
lon2 <- lon[n <= nlon/2];
lon_transf <- c(lon1, lon2);

tas_transf <- array(dim=c(nlon, nlat, Nt));
tas_transf[n <= nlon/2,,] <- tas[n > nlon/2,,];
tas_transf[n > nlon/2,,] <- tas[n <= nlon/2,,];

library(rasterVis);

tasMap <- raster(ncol=320, nrow=160, xmn=min(lon_transf), xmx=max(lon_transf), ymn=min(lat), ymx=max(lat));
tasMap2 <- raster(ncol=320, nrow=160, xmn=min(lon_transf), xmx=max(lon_transf), ymn=min(lat), ymx=max(lat));
tasMap3 <- raster(ncol=320, nrow=160, xmn=min(lon_transf), xmx=max(lon_transf), ymn=min(lat), ymx=max(lat));
tasMap4 <- raster(ncol=320, nrow=160, xmn=min(lon_transf), xmx=max(lon_transf), ymn=min(lat), ymx=max(lat));

for(i in seq(1,nlat)){
 
  tasMap[(nlon*(i-1)+1):(i*nlon)] <- tas_transf[seq(1, nlon),(nlat-i)+1,1];
 
}

for(i in seq(1,nlat)){
 
  tasMap2[(nlon*(i-1)+1):(i*nlon)] <- tas_transf[seq(1, nlon),(nlat-i)+1,10];
 
}

for(i in seq(1,nlat)){
 
  tasMap3[(nlon*(i-1)+1):(i*nlon)] <- tas_transf[seq(1, nlon),(nlat-i)+1,100];
 
}

for(i in seq(1,nlat)){
 
  tasMap4[(nlon*(i-1)+1):(i*nlon)] <- tas_transf[seq(1, nlon),(nlat-i)+1,500];
 
}

tasMapTot <- stack(tasMap, tasMap2, tasMap3, tasMap4);
names(tasMapTot) <- c("RCP2.6", "RCP4.5", "RCP6.0", "RCP8.5");


library(grDevices);
library(maps);

couleur <- c(rgb(255/255,255/255,178/255), rgb(254/255,204/255,92/255), rgb(253/255,141/255,60/255), rgb(227/255,26/255,28/255));
color_plate <- colorRampPalette(couleur, space="rgb");

brks <- seq(min(tas), max(tas), length.out=100);
Ncol <- length(brks);

options(device=X11());
spplot(tasMapTot, col.regions=color_plate(Ncol-1), scales=list(draw=T));

savePlot(filename="essai_raster.png", type="png");



J'ai essayé d'ajouter les continents avec les fonctions map(add=T) et addland() mais j'obtiens l'erreur suivante:

> map(add=T);
Error in plot.xy(xy.coords(x, y), type = type, ...) :
plot.new has not been called yet

Quelqu'un aurait une idée?

Encore merci pour l'aide :-)

[Renaud Lancelot]

Le package raster est vraiment très simple à utiliser (et peut gérer des rasters au format netcdf). Ce n'est pas un package réservé au traitement des images satellites même s'il facilite beaucoup les choses.

Votre code pourrait être simplifié par l'utilisation de ce package. Noter aussi que les points virgules en fin de ligne ne servent à rien. Voici ce que je ferais en pratique (ça marche avec les données que vous avez fournies):

Code : Tout sélectionner

library(raster)
## Crée des rasters
for(i in seq(4)){
    assign(paste0("r", i),
           flip(raster(x = t(tas_transf[ , , i + 10 * (i - 1)]),
                       xmn = -160, xmx = 160, ymn = -90, ymx = 90),
                direction = "y"))
    }
## Empile les rasters
r <- stack(r1, r2, r3, r4)
## Nomme les rasters (à adapter)
names(r) <- paste("Date", 1:4)

## Cartes "simple"
spplot(r,
       col.regions = color_plate(100))

## Cartes avec pays
library(maps)
world <- map(regions = ".", fill = TRUE, plot = FALSE)
library(maptools)
## Conversion map --> SpatialPolygons
nam <- strsplit(world$names, ":")
id  <- sapply(nam, function(x) x[1])
World <- map2SpatialPolygons(world, IDs = id, proj4string = CRS("+proj=longlat"))

## Trace les cartes
spplot(r,
       col.regions = color_plate(100),
       panel = function(...){
           panel.gridplot(...)
           sp.polygons(World)
           })


[Renaud Lancelot]

J'ai essayé d'ajouter les continents avec les fonctions map(add=T) et addland() mais j'obtiens l'erreur suivante:

> map(add=T);
Error in plot.xy(xy.coords(x, y), type = type, ...) :
plot.new has not been called yet

Quelqu'un aurait une idée?

Encore merci pour l'aide :-)

Voir ma réponse précédente. L'idée est de transformer les objets de classe "map" en objets spatiaux de classe compatible avec celle produite par raster: fonction map2SpatialPolygons du package maptools.

[Quynh-Mai Chu]

Merci beaucoup Renaud, ce code fut très instructif. Je ne connaissais pas les fonctions assign et flip (qui raccourcissent beaucoup de lignes de codes) :-).