[R-sig-Geo] Execute Extract function on several Raster with points layer
Bede-Fazekas Ákos
b|@|ev||@t @end|ng |rom gm@||@com
Thu Sep 17 19:57:18 CEST 2020
Hello Gaëtan,
so as far as I understand, you have 3 main folders:
"Max_T", ? and ?
and in alll the three folders, there are subfolders
"1961", "1962", ... "1970"
In each folder, there are 366 raster files, for which the file naming
conventions are not known by us, but some of the files are called
"max1961_1.asc", "max1961_2.asc", ... "max1961_366.asc" (in case of
T_max and year 1961)
In this case, the 10980 layer that belongs to T_max can be read to one
large RasterStack in this way:
tmax_filenames <- c(outer(X = as.character(1:366), Y =
as.character(1961:1970), FUN = function(doy, year) paste0("N:/400074
Conservation des sols et CC/Data/Climate data/Climate-10km/Max_T/",
year, "/max", year, "_", doy, ".asc")))
tmax_raster <- stack(tmax_filenames)
You can give self-explanatory names to the raster layers:
names(tmax_raster) <- c(outer(X = as.character(1:366), Y =
as.character(1961:1970), FUN = function(doy, year) paste0(year, "_", doy)))
But if the structure of the rasters are the same (i.e. the cell size,
extent, projection), then I recommend you to do the raster-vector
overlay once, save the cell numbers that you are interested in, and then
in nested for loops (one loop for the climate variable, one for the year
and one for the day) read the rasters one-by-one, extract the values
according to the cell numbers, and save the result in a previously
created data.frame. In this way, you may not encounter memory issues.
Although, it will take a lot of time...
HTH,
Ákos Bede-Fazekas
Hungarian Academy of Sciences
2020.09.17. 19:28 keltezéssel, Gaetan Martinelli írta:
> Hello everyone, R team,
>
> Sorry in advance for this long message. Your help will be invaluable.
>
> For a few days now i have been blocked to execute a task on R. I will
> try to synthesize my problem.
>
> I have several raster. I have an ASCII file for each day of a year
> with a single band. For 30 years, and for three climatic variables on
> grid 10km/10km (T_min, T_max, Precipitation). So i have a total around
> of 32 940 raster files (366days*30years*3variables).
>
> Also, i have a layer of aroud 1000 points.
>
> I tried to use the Stack function and then make the intersection for
> each raster files with my 1000 points.
> I cannot create an independent matrix for all my files where i applied
> the "extract" function, to then concatenate all my matrices in order
> to have a single table.
>
> I tried this, exemple for 10 years et only T_Max (my files are
> organized the same for my two other variables) :
> *#Datapoints*
> Datapoints<-readOGR(dsn="H:/Inventaire/R/final",
> layer="Centroid_champs")
> Datapoints<- spTransform (Datapoints, CRS ("+init=epsg:4326") ) # 1022
> points in the data
> st_crs(Datapoints)
> *#Rasters files*
> folders = list(
> file.path('N:','Data','Climate data','Climate-10km','Max_T','1961'),
> #Each year includes daily data, the names of my several raster is
> "max1961_1", "max1961_10", "max1961_100", etc...
> file.path('N:','Data','Climate data','Climate-10km','Max_T','1962'),
> file.path('N:','Data','Climate data','Climate-10km','Max_T','1963'),
> file.path('N:','Data','Climate data','Climate-10km','Max_T','1964'),
> file.path('N:','Data','Climate data','Climate-10km','Max_T','1965'),
> file.path('N:','Data','Climate data','Climate-10km','Max_T','1966'),
> file.path('N:','Data','Climate data','Climate-10km','Max_T','1967'),
> file.path('N:','Data','Climate data','Climate-10km','Max_T','1968'),
> file.path('N:','Data','Climate data','Climate-10km','Max_T','1969'),
> file.path('N:','Data','Climate data','Climate-10km','Max_T','1970')
> )
> files = unlist(sapply(folders, function(folder) {
> list.files(folder, full.names=TRUE)
> }))
> files
>
> MET <- lapply(files, raster)
> s <- raster::stack(MET)
>
> output <- list()
> for(i in 1:length(MET)){
> output[[i]] <- extract(s, Datapoints)
> names(output)[[i]] <- paste("Année", MET[i], sep = "_")
> }
> Also, i tried that :
> p1 <- 1022 (ID of my DataPoints) ; p2 <- 1 (column where there are the
> values extracted from my raster) ; p3 <- 3660 # 3660matrix (366
> day* 10 years)
> matlist <- list(array(NA,c(p1,p2,p3))) # doing a list of independant
> matrix
>
> for(i in seq_along(MET)){
>
> matlist[[i]] <- extract(s, Datapoints)
> }
> But, nothing works...
> I would like my script to perform these actions :
> - For each Raster in my Rasterstack, extract the climatic data values
> and link them to my "Datapoints",
> - Take the name of my file, take the first three characters of the
> name to get a column of my weather variable, here, "T_Max" (column
> with my raster values) ; Take the following four characters then
> report this information in a new column "Year", and finally, take the
> last characters of the file name to create a new column "Day".
> - Concatenate all the independent output matrices corresponding to
> each intersection made with my different raster files
> In the end, I would have a huge table, but one that will allow me to
> do my analysis :
> Table with 9 attributes (6 attributs of my points + Year + Day +
> T_Max) like this :
> ID Datapoint Year Day T_Max
> 1 1960 1
> 2 1960 1
> …... 1960 1
> 1022 1960 1
> 1 1960 2
> 2 1960 2
> …... 1960 2
> 1022 1960 2
> ….. ….. …..
> 1 1970 1
> 2 1970 1
> …... 1970 1
> 1022 1970 1
> 1 1970 2
> 2 1970 2
> …... 1970 2
> 1022 1970 2
> ….. ….. …..
>
> Could a loop do this task ?
>
> I'm sorry, i am gradually learning to manipulate R, but this exercise
> is more difficult than expected...Please feel free to tell me if my
> question is inappropriate.
>
> Thank you very much in advance for your answers. Your help or your
> comments will be really appreciated.
>
> Have a good day.
>
> Gaëtan Martinelli
> Water and Agriculture research professional in Quebec.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list