Extraction de valeurs le long de bandes circulaires personnalisées (distance radiale constante) à partir d'un raster

Question

J'ai un grand objet raster importé à partir d'un fichier image .tif. Je veux analyser les pixels de ce raster, à une distance radiale donnée du centre, afin d'identifier certains phénomènes axisymétriques que je peux remarquer dans l'image (ci-dessous). Pour ce faire, je veux extraire les valeurs de pixels qui se croisent avec des bandes circulaires d'un rayon donné (et de la largeur) à partir du centre de l'image (indiqué dans la figure ci-dessous).

J'ai exploré plusieurs options pour ce faire, y compris la fonction d'extraction et le paquet d'imageur. Dans le paquetage imageur, vous pouvez facilement extraire des valeurs le long de lignes ou de colonnes, mais je n'ai pas trouvé de fonctions permettant d'extraire des valeurs avec des formes personnalisées comme je le demande. Je pourrais utiliser la fonction d'extraction avec l'argument SpatialPolygons, mais pour cela, je devrai entrer les emplacements de tous les points dans un objet polygone, ce qui nécessiterait une localisation très élevée des points (depuis les points, je pense , sont joints par des éléments de ligne). De plus, je veux faire varier la densité numérique des bandes sur lesquelles j'extrais les valeurs des pixels (et plus tard la moyenne) donc cette méthode serait à la fois fastidieuse et inflexible. Je me demandais donc si certains d'entre vous avaient des suggestions pour s'attaquer à ce problème.

> str(imRaster) 
Formal class 'RasterLayer' [package "raster"] with 12 slots 
    ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots 
    .. .. ..@ name  : chr "C:\\Users\\Nandu\\input_images\\All\\101_1A_1000ms.tif" 
    .. .. ..@ datanotation: chr "INT2U" 
    .. .. ..@ byteorder : chr "little" 
    .. .. ..@ nodatavalue : num -Inf 
    .. .. ..@ NAchanged : logi FALSE 
    .. .. ..@ nbands  : int 1 
    .. .. ..@ blockrows : int 1 
    .. .. ..@ blockcols : int 1392 
    .. .. ..@ driver  : chr "gdal" 
    .. .. ..@ open  : logi FALSE 
    ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots 
     . 
     . 
     . 

Alternativement, s'il existe d'autres méthodes (autres que par l'importation d'une image .tif à raster) qui permettrait une telle opération, qui pourrait aussi être utile. S'il vous plaît, faites-moi savoir.

Merci d'avance!

EDIT: Afin de rendre le problème reproductible, j'ai ajouté ici le lien vers l'image que je veux utiliser. Il est une image tiff, et peut être téléchargé ici

Link to the tiff image that I am trying to process in R

library(tiff) 
library(raster) 
imRaster = raster(file_path to the image) 
plot(imRaster) 
xy = cbind(684.4228, 599.0458) # Gives a rough location of the center in the image coordinates 

Ce code peut être exécuté afin d'obtenir la trame décrite ci-dessus et le visualiser. J'espère que cela aide!

Answer 1

Commençons par le charger la trame et l'affectation d'un système de référence de coordonnées faux (CRS). Nous utilisons un métrique de sorte que 1 pixel dans votre image se réfère à un "mètre" dans le pseudo pseudo raster.

library(raster) 

imRaster <- raster("~/Downloads/46-ECN300448-216.tif") 
imRaster@crs <- CRS("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs") 

plot(imRaster, col = grey.colors(255)) 

Ensuite, comme vous l'avez mentionné, nous utilisons votre point central xy et le convertir en un SpatialPoint

xy <- cbind(684.4228, 599.0458) 
xySp <- SpatialPoints(coords = xy) 
plot(xySp, add = T) 

Ensuite, nous utilisons le paquet rgeos pour créer deux tampons radiaux autour du point . Dans cet exemple, nous utilisons 200 et 300 px/m.

library(rgeos) 

outer <- 300 
inner <- 200 

bufOut <- gBuffer(xySp, width = outer) 
bufIn <- gBuffer(xySp, width = inner) 

strip <- bufOut - bufIn 
plot(strip, add = T, col = "#FF000050") 

Enfin, nous pouvons utiliser le strip pour masquer l'image raster et calculer des statistiques (ou même utiliser les valeurs d'origine) avec getValues().

m <- mask(imRaster, mask = strip) 

# plot(m) # plot the mask if you want to see what it looks like 

mean(getValues(m), na.rm = T) 
# [1] 17004.24 

Remarque importante: Seules les cellules ayant leur centre dans le polygone sont restent après masquage. Vous pouvez aborder ce problème en jouant avec les tampons.