Suite

Interpolation spatio-temporelle dans R ou ArcGIS ?


J'essaie de calculer la valeur moyenne des précipitations à partir d'un certain nombre de points à l'aide de l'outil Distance pondérée inverse dans ArcGIS 9.3.

Mon problème est que : chaque point a sa propre série temporelle, donc le processus d'interpolation devrait pouvoir s'effectuer pour toutes les années (sorte d'itération pour ainsi dire).

Voici un exemple de table attributaire :

ID X Y Nom Rain1990 Rain1991 Rain1992 Rain1993… Rain2010 1 xx1 yy1 AA 1210 1189 1863 1269… 2 xx2 yy2 BB 1492 1502 2187 1923…

Quelqu'un pourrait-il me montrer comment faire?


Edit 1 : j'ai finalement fait cela en utilisant du code C++ qui nécessitait une grille de masque ArcGIS, des fichiers de données et des emplacements de tous les points.


Edit 2: J'ai récemment utilisé R pour faire cette tâche d'interpolation. Vous pouvez utiliser soithydroTSM,gstatou alorsespace-tempspaquets. Quelques exemples de liens ci-dessous :

http://spatial-analyst.net/wiki/index.php?title=Spatial_interpolation_exercises_%28NL%29

http://www.geostat-course.org/Topic_Bivand_2012


Edit 3: Ajout d'un exemple de travail ci-dessous pour les futurs lecteurs


J'ai résolu ce problème en insérant un itérateur "Feature Selection" dans un modèle. (Dans la fenêtre ModelBuilder, sous le menu Insertion->Itérateurs.)

Utilisez votre champ de temps comme variable de « groupe par ». En faisant cela, le modèle itérera une fois pour chaque fois dans votre classe d'entités.

Ensuite, attachez votre outil d'interpolation préféré (spline, IDW, peu importe) à la sortie de l'entité de l'itérateur. Exécutez le modèle, partez en vacances pendant quelques semaines, et à votre retour, vous aurez autant de grilles que vous avez de points temporels dans la classe d'entités.

Notez que cette solution suppose que vous avez des points d'échantillonnage temporels discrets avec une date ou un champ numérique qui indique un point temporel unique pour chaque enregistrement de votre ensemble d'entités. Si vous utilisez le format "heure de début" et "heure de fin", ce n'est peut-être pas aussi simple.


Il semble que ce fil soit répondu par l'outil IDW, mais si vous deviez demander et saisir l'année de début, puis parcourir les champs de l'année à l'aide d'une variable en ligne dans le générateur de modèles, ce serait un moyen plus élégant de gérer la modélisation .

PS: je suis d'accord avec @AndyW que si vous le résolvez en utilisant l'IDW, postez vous-même comme réponse, puis "cochez avec la coche"


Ajouter ma propre solution en utilisantR& données de précipitations aléatoires

library(tidyverse) library(sp) # pour les coordonnées, CRS, proj4string, etc. library(gstat) library(maptools) # Coordonnées des cellules de précipitation maillées precGridPts <- ("ID lat long 1 46.78125 -121.46875 2 46.84375 -121.53125 3 46.84375 - 121.46875 4 46.84375 -121.40625 5 46.84375 -121.34375 6 46.90625 -121.53125 7 46.90625 -121.46875 8 46.90625 -121.40625 9 46.90625 -121.34375 10 46.90625 -121.28125 11 46.96875 -121.46875 12 46.96875 -121.40625-12 46.96875 46.96875 -121.15625 ") # Lire les cellules de précipitation precGridPtsdf <- read.table(text = precGridPts, header = TRUE)

Convertir en objet sp

sp::coordinates(precGridPtsdf) <- ~long + lat # longitude en premier

Ajoutez un système de référence spatiale (SRS) ou un système de coordonnées de référence (CRS).

# Base de données CRS : http://spatialreference.org/ref/epsg/ sp::proj4string(precGridPtsdf) <- sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84") str(precGridPtsdf) # > Classe formelle 'SpatialPointsDataFrame' [package "sp"] avec 5 slots #>… @ data :'data.frame' : 16 obs. de 1 variable : #>… $ ID : int [1:16] 1 2 3 4 5 6 7 8 9 10… #>… @ coords.nrs : int [1:2] 3 2 #>… @ coords : num [1:16, 1:2] -121 -122 -121 -121 -121… #>… - attr(*, "dimnames")=Liste de 2 #>… $ : chr [1:16] "1" "2" "3" "4"… #>… $ : chr [1:2] "long" "lat" #>… @ bbox : num [1:2, 1:2] -121,5 46,8 -121,2 47 # >… - attr(*, "dimnames")=Liste de 2 #>… $ : chr [1:2] "long" "lat" #>… $ : chr [1:2] "min" "max" # >… @ proj4string:Classe formelle 'CRS' [package "sp"] avec 1 slot #>… @ projargs: chr "+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"

Convertir en UTM 10N

utm10n <- "+proj=utm +zone=10 ellps=WGS84" precGridPtsdf_UTM <- spTransform(precGridPtsdf, CRS(utm10n))

Données hypothétiques sur les précipitations annuelles générées à l'aide de la distribution de Poisson.

precDataTxt <- ("ID PRCP2016 PRCP2017 PRCP2018 1 2125 2099 2203 2 2075 2160 2119 3 2170 2153 2180 4 2130 2118 2153 5 2170 2083 2179 6 2109 2008 2107 7 2109 2189 2093 8 2058 2170 2067 9 2154 2119 2139 10 2056 2184 2120 11 2080 2123 2107 12 2110 2150 2175 13 2176 2105 2126 14 2088 2057 2199 15 2032 2029 2100 16 2133 2108 2006" ) precData <- read_table2(precDataTxt, col_types = cols(ID = "i"))

Fusionner le cadre de données Prec avec le fichier de formes Prec

precGridPtsdf <- merge(precGridPtsdf, precData, by.x = "ID", by.y = "ID") precdf <- data.frame(precGridPtsdf)

Fusionner le cadre de données des précipitations avec le fichier de formes des précipitations (UTM)

precGridPtsdf_UTM <- merge(precGridPtsdf_UTM, precData, by.x = "ID", by.y = "ID") # échantillon étendue region_extent <- structure(c(612566.169007975, 5185395.70942594, 639349.654465079, 5205871.0782451),(2.D , 2L), .Dimnames = list(c("x", "y" ), c("min", "max")))

Définissez l'étendue de l'interpolation spatiale. Agrandissez-le de 4 km dans chaque direction

x.range <- c(region_extent[1] - 4000, region_extent[3] + 4000) y.range <- c(region_extent[2] - 4000, region_extent[4] + 4000)

Créez la grille souhaitée à une résolution de 1 km

grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 1000), y = seq(from = y.range[1], to = y .range[2], by = 1000)) # Convertir la grille en coordonnées d'objets spatiaux (grd) <- ~x + y # Utiliser la même projection que bound_UTM proj4string(grd) <- "+proj=utm +zone=10 ellps =WGS84 +ellps=WGS84" gridded(grd) <- TRUE

Interpoler à l'aide de la pondération inverse de la distance (IDW)

idw <- idw(formula = PRCP2016 ~ 1, locations = precGridPtsdf_UTM, newdata = grd) #> [interpolation pondérée par la distance inverse] # Nettoyer idw.output = as.data.frame(idw) names(idw.output)[1 :3] <- c("Longitude", "Latitude", "Précipitations") precdf_UTM <- data.frame(precGridPtsdf_UTM)

Tracer les résultats de l'interpolation

idwPlt1 <- ggplot() + geom_tile(data = idw.output, aes(x = Longitude, y = Latitude, fill = Precipitation)) + geom_point(data = precdf_UTM, aes(x = long, y = lat, size = PRCP2016 ), shape = 21, color = "red") + viridis::scale_fill_viridis() + scale_size_continuous(name = "") + theme_bw() + scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(expand = c( 0, 0)) + theme(axis.text.y = element_text(angle = 90)) + theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0))) idwPlt1

### Maintenant en boucle chaque année list.idw <- colnames(precData)[-1] %>% set_names() %>% map(., ~ idw(as.formula(paste(.x, "~ 1" )), locations = precGridPtsdf_UTM, newdata = grd)) #> [interpolation pondérée en distance inverse] #> [interpolation pondérée en distance inverse] #> [interpolation pondérée en distance inverse] idw.output.df = as.data.frame(list. idw) %>% as.tibble() idw.output.df #> # Un tibble : 1 015 x 12 #> PRCP2016.x PRCP2016.y PRCP2016.var1.pred PRCP2016.var1.var PRCP2017.x #> *      #> 1 608566. 5181396. 2114. NA 608566. #> 2 609566. 5181396. 2115. NA 609566. #> 3 610566. 5181396. 2116. NA 610566. #> 4 611566. 5181396. 2117. NA 611566. #> 5 612566. 5181396. 2119. NA 612566. #> 6 613566. 5181396. 2121. NA 613566. #> 7 614566. 5181396. 2123. NA 614566. #> 8 615566. 5181396. 2124. NA 615566. #> 9 616566 . 5181396. 2125. NA 616566. #> 10 617566. 5181396. 2125. NA 617566. #> #… avec 1005 lignes supplémentaires et 7 variables supplémentaires : PRCP2017.y , #> # PRCP2017.var1.pred , PRCP2017.var1.var , PRCP2018.x , #> # PRCP2018.y , PRCP2018.var1.pred , PRCP2018.var1.var 


Voir la vidéo: Prediction of LULC changes for species distribution modeling using cellular automata and ANN (Octobre 2021).