1

I am trying to select pixels from a MODIS raster (with a 250 meters resolution and covering all Europe) based on their location. For that I am using a shapefile that contains ca. 200 points (also distributed across all Europe) representing the location of different time series. I am particularly interested in selecting all pixels that intersect with a 1000m-buffer around each point. Finally, for further analysis, I need that the selected pixels (those that intersect with the buffers) will have the id of the point (i.e., the id of the time series). To save time I am trying to avoid working with vector data. So I have rasterized the shapefile using the extension and resolution of the MODIS raster and ratify it as suggested here. The problem I have is that after using raster::buffer the RAT (raster attribute table) disappears and all pixel has a value = 1. Does anyone know how can I buffer a raster and keep the RAT?

Below a reproducible code showing the problem. The layers are subset form my data and are accessible here:

library(raster) library(rgdal) library(rasterVis) # load Alps Modis MODIS_Alps <- raster("/ModisAlpsCrop.tif") plot(MODIS_Alps) # reproject sr <- "+proj=longlat +datum=WGS84" MODIS_Alps_LatLong <- projectRaster(MODIS_Alps, crs = sr) plot(MODIS_Alps_LatLong) # load TS BeechAlps <- read_csv("/BeechAlps.csv") BeechAlps <- BeechAlps %>% st_as_sf(coords = c("X", "Y"), crs = 4326) %>% dplyr::select(ID) plot(BeechAlps, col= "black", add=TRUE) # get names nam <- unique(BeechAlps$ID) # create a data.frame nam_df <- data.frame(ID2 = 1:length(nam), nam = nam) # Place IDs BeechAlps$ID2 <- nam_df$ID2[match(BeechAlps$ID,nam_df$nam)] # rasterize MODIS_Alps_LatLong r <- raster(ncol = 185, nrow = 162, ext = extent(MODIS_Alps_LatLong)) # resolution and extent from MODIS ras <- rasterize(x = BeechAlps, y = r, field = "ID2") # ratify raster r2 <- ratify(ras) # Create levels rat <- levels(r2)[[1]] rat$names <- nam_df$nam rat$IDs <- nam_df$ID2 rat levels(r2) <- rat rasterVis::levelplot(r2) r2 # buffer b <- buffer(r2, width=1000) b plot(b, col="black")
Improve this question

2 Answers 2

Reset to default
0

Am I missing something evident or is this whole concept of raster attribute tables only existing inside ESRI environment? I always thought that only features can be characterized by further attributes, but raster data is simply limited to the storage of one value per grid cell. Hence, raster attribute tables sound like something artificial trying to combine the best of both worlds using maybe some kind of a join on the raster values to be able to append further information.

However, and excuse me if my approach does not suit you because I'm not respecting these boundary conditions. Nevertheless, I would like to provide an alternative solution, which seems more straightforward, at least from my point of view:

library(readr) library(dplyr) library(terra) #> terra 1.6.3 library(sf) #> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE # Import MODIS data as SpatRaster r <- rast("ModisAlpsCrop.tif") # Set coordinate reference systems for next step crs_4326 <- st_crs(4326) crs_laea <- crs(r) # Create data frame from csv file and tidy column names for join later on data <- read_csv("BeechAlps.csv") |> mutate(name = ID, ID = row_number()) # Transform df to sf in WGS 84, reproject to LAEA, transform to SpatVector p <- st_as_sf(data, coords = c("X", "Y"), crs = crs_4326) |> st_transform(crs_laea) |> vect() # Construct buffer polygons around your points b <- buffer(p, width = 1000) # Inspect progress so far plot(r) plot(p, add = TRUE) plot(b, add = TRUE) 

 # Extract values from raster using buffer polygons without aggregation via `fun = `. # Benchmarking via `tictoc` reveals this operation was carried out in 0.11 s on my machine. results <- extract(r, b, na.rm = TRUE) # Join station names to the results based on ID column results <- inner_join(x = results, y = data, by = "ID") |> select(ID, name, ModisAlpsCrop) head(results) #> ID name ModisAlpsCrop #> 1 1 FASYCH0038A 1041.3146 #> 2 1 FASYCH0038A 1628.9230 #> 3 1 FASYCH0038A 1829.2584 #> 4 1 FASYCH0038A 2088.8303 #> 5 1 FASYCH0038A 725.1849 #> 6 1 FASYCH0038A 1101.9822

Some thoughts and further explanations:

  • Since rgdal will retire by the end of 2023, I switched from raster to its successor terra.

  • Unless performance is not a restricting issue in real terms, I would not transform vector to raster and vice versa for the purpose of optimization, because you simply run into several limitations - as you did - and also your code becomes more complex making it less maintainable over time. You can benchmark once making use of e.g. tictoc or with n iterations as a batch job for statistical purposes using microbenchmark.

  • It's always a good idea to use common coordinate reference systems before performing any kind of operation on spatial data. But it's a better idea to reproject your vector data instead of your raster data to prevent information blurring due to resampling. Also, a projected reference system makes it easier to construct buffer polygons with a specific width in meters.

  • Your main problem, as far as I got it, is that you lost information from your rasterized buffers so that you ran into issues because there was no relation evident between the extracted grid cells and the station ID used. Since order is strictly maintained when using extract, names stored in data can be joined on the results by ID, c.f. results. Clean and easy from my point of view.

So on the one hand, this can be regarded as a more or less minimalistic solution to your problem. On the other hand, I was not really able to answer your initial question considering RAT and used a completely different approach. Feel free to decide if this suits you nevertheless.

Improve this answer
2
  • 1
    RATs are certainly not an Esri concept. e.g. the KEA format, based on HDF5, has excellent support for RATs. github.com/ubarsc/kealib Commented Aug 4, 2022 at 0:22
  • Thanks for the heads-up! Did not have much contact with RATs yet - as you may have noticed - but will look into this. Commented Aug 4, 2022 at 1:29
0

@falk-env, your code has helped me a lot. I was missing the coordinates for each MODIS pixel within the buffer. I have just added a few lines to your code and get what I needed.

library(readr) library(dplyr) library(terra) #> terra 1.6.3 library(sf) #> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE # setwd("G:/My Drive/RAT") # Import MODIS data as SpatRaster r <- rast("ModisAlpsCrop.tif") # Set coordinate reference systems for next step crs_4326 <- st_crs(4326) crs_laea <- crs(r) # Create data frame from csv file and tidy column names for join later on data <- read_csv("BeechAlps.csv") %>% mutate(name = ID, ID = row_number()) # Transform df to sf in WGS 84, reproject to LAEA, transform to SpatVector p <- st_as_sf(data, coords = c("X", "Y"), crs = crs_4326) %>% st_transform(crs_laea) %>% vect() p # Construct buffer polygons around your points b <- buffer(p, width = 1000) # rasterize buffer r r2 <- rast(ncol = 144, nrow = 172, # resolution from MODIS xmin = 4178329, xmax = 4218405, ymin = 2676126, ymax = 2709246) # extent from MODIS ras <- rasterize(x = b, y = r2, field = "ID") plot(ras) ras # Inspect progress so far plot(r) plot(ras, add = TRUE) plot(b, add = TRUE) # convert rasterized buffer into a data frame ras_df <- as.data.frame(ras, xy=TRUE) summary(ras_df) # sort the data frame to later add the IDs #(extract function will only work with two columns - X Y) ras_df <- ras_df[order(ras_df$ID),] summary(as.factor(ras_df$ID)) ID <- ras_df$ID # remove ID ras_XY <- ras_df[,-3] # extract results <- extract(r, ras_XY, xy= TRUE, na.rm = TRUE) head(results) nrow(results) == nrow(ras_XY) nrow(results) == length(ID) results <- cbind(results, ID) # Join station names to the results based on ID column results <- left_join(x = results, y = data, by = "ID") %>% select(ID, name, x, y) head(results)

Thank you once again!

Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.