Suite

Est-il possible de découper un fichier de formes dans une image () dans R ?


Je travaille sur la cartographie de certaines données netcdf dans R et j'ai eu du mal à essayer d'effectuer un clip avec un fichier de formes de pays importé que j'utilise. Voici ma question d'hier. Si vous savez quelque chose à ce sujet, n'hésitez pas à venir et à répondre à cette question. :)

MAIS je vais recommencer ici. Donc, pour mon travail, nous prévoyons de cartographier des tonnes et des tonnes de données météorologiques à travers le monde. Pour s'entraîner, le gars qui a fourni nos données nous a fourni des "exemples" de données avec lesquelles jouer. Mon patron aimerait que je crée mes cartes en R. Ces cartes sont destinées à être des cartes de période de retour. Les données sont au format netcdf. J'ai pu mapper les données netcdf avec la fonction levelplot du package rasterVis. J'ai également pu mapper les données netcdf avec une simple fonction d'image. Voici mon code (désolé pour les conventions de nommage funky ; je viens de jouer avec ça):

bibliothèque(maptools) bibliothèque(ncdf4) bibliothèque(RCorBrewer) bibliothèque(raster) setwd("D:stuff") #Madagascar shapefile madagascar <- readShapeSpatial("MDG_adm0.shp") plot(madagascar) #Appeler sur le fichier netcdf ncdf .data <- nc_open("swio_rpmaps_200_83.nc") #Appel sur les périodes de retour du vent de 10 ans retourne <- ncvar_get(ncdf.data, "y010") #Définir la latitude et la longitude lon <- ncdf.data$dim$longitude$ vals lat <- ncdf.data$dim$latitude$vals #plot image of windspeeds image(lon, lat, return, col = cm.colors(9,alpha=.6), add = TRUE) #plot Madagascar on top of graphique des vitesses du vent (madagascar, ajouter = T)

Voici l'image du résultat :

Existe-t-il un moyen de découper le fichier de formes Madagascar sur cette image ? Ou une autre façon de manipuler les données pour que je puisse réaliser un clip ? Idéalement, je n'aurais que les données associées dans les limites du pays que je cartographie.


Reprenant la suggestion de @hrbrmstr voici le code. Puisque vous ne fournissez pas d'exemple de fichier, j'essaie d'imiter votre situation et de créer un exemple de la manière dont vous pourriez le faire avec votre fichier netcdf. Le shapefile que j'ai extrait de Natural Earth.

library(maptools) library(raster) madagascar <- readShapeSpatial("MDG_adm0.shp") ## make up some values ​​# J'essaie de reproduire la façon dont vous le feriez avec le fichier netcdf… lon <- seq(42, 51 , .05) lat <- seq(-29, -10, .1) renvoie <- runif(length(lat) * length(lon)) image(lon, lat, matrix(returns, length(lon), length( lat)), col = cm.colors(9,alpha=.6)) plot(madagascar, add = T)

Ensuite, nous utilisons les valeurs du fichier netcdf pour créer un objet raster réel (un RasterLayer pour être exact) avec lequel vous pouvez travailler. Pour ce faire, nous devons définir les dimensions du raster (nrows et ncols) ainsi que les valeurs min et max pour les deux dimensions (x et y). On utilise alors lesetValuesfonction pour ajouter les valeurs (dans votre cas, les vitesses du vent) pour chaque cellule raster.

METTRE À JOUR: Au lieu de faire cela, créez un objet Raster directement comme indiqué par @RobertH .

# créer l'objet raster en utilisant les valeurs ci-dessus r <- raster(nrows=length(lon), ncols=length(lat), xmn=min(lon), xmx=max(lon), ymn=min(lat) , ymx=max(lat)) r <- setValues(r, return) plot(r) plot(madagascar, add = T)

Enfin, nous utilisons lemasquerfonction pour découper le raster. Cette fonction prend n'importe quel objet Spatial* comme masque, un SpatialPolygonsDataframe (comme celui de Madagascar) étant l'un d'entre eux.

#clip r.clipped <- mask(r, madagascar) plot(r.clipped)


Ajout à la réponse de cengel, pour montrer comment accéder au fichier ncdf en tant qu'objet RasterBrick, puis continuer :

library(raster) madagascar <- shapefile("MDG_adm0.shp") b <- raster("swio_rpmaps_200_83.nc", var= "y010") bb <- crop(b, madagascar) bb <- mask(bb, madagascar)