0

I'm uncertain as to why my geometry coordinates are empty.

The dataset I have derived the X and Y coordinates had the same CRS projection as the MODIS data landcover. The following data is MODIS MCD12Q1 yearly landcover data taken from the extent of 17 BCR regions in North America.

Given as:

landcover
class      : RasterStack 
dimensions : 7149, 10258, 73334442, 10  (nrow, ncol, ncell, nlayers)
resolution : 463.3127, 463.3127  (x, y)
extent     : -9793968, -5041306, 2716866, 6029088  (xmin, xmax, ymin, ymax)
crs        : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
names      : y2010, y2011, y2012, y2013, y2014, y2015, y2016, y2017, y2018, y2019 
min values :     0,     0,     0,     0,     0,     0,     0,     0,     0,     0 
max values :   255,   255,   255,   255,   255,   255,   255,   255,   255,   255 

Here is my dataset with the following X and Y coordinates:

 A tibble: 6 x 8
     X1 encounter    id layer  water  year         X        Y
  <dbl>     <dbl> <dbl> <dbl>  <dbl> <dbl>     <dbl>    <dbl>
1     1    0.0350     1     1 0.0714  2011  -7158197 5814408
2     2    0.0507     2     1 0.0714  2011  -7158197 5812328
3     3    0.0310     3     1 0.111   2011  -7657517 5793608
4     4    0.0585     4     1 0.111   2011  -7661897 5791528
5     5    0.0437     5     1 0.111   2011  -7659707 5791528
6     6    0.0590     6     1 0.125   2011  -7668467 5783208

I proceeded with the following code to reshape the geometry points and corresponding values in the column into a long dataframe format:

bird_buff <- ebird_buff$data[1] %>%
  bind_cols %>%
  st_cast(to = "POINT") %>%
  dplyr::mutate(
    X =  sf::st_coordinates(geometry)[,1], #retrieve X coord
    Y =  sf::st_coordinates(geometry)[,2]  #retrieve Y coord
  ) %>%
  sf::st_drop_geometry()


bird_buff$year <- as.integer(bird_buff$year)
bird_buff$X <- as.numeric(bird_buff$X)
bird_buff$Y <- as.numeric(bird_buff$Y)

agg_factor <- round(2 * neighborhood_radius / res(landcover))
r <- raster(landcover) 
bird.raster <- bird_buff %>% st_as_sf(coords = c("X", "Y")) %>% rasterize(r) 

output:

class      : RasterBrick 
dimensions : 1430, 2052, 2934360, 3  (nrow, ncol, ncell, nlayers)
resolution : 2316.564, 2316.564  (x, y)
extent     : -9793968, -5040379, 2716402, 6029088  (xmin, xmax, ymin, ymax)
crs        : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
source     : memory
names      : ID, year, encounter 
min values : NA,   NA,        NA 
max values : NA,   NA,        NA 

To download the files:

file1 <- "modis_mcd12q1_umd_2010.tif"
file2 <- "shape.shp"
file3 <- "shape.shx"


dir.create("data", showWarnings = FALSE)
if (!file.exists(file.path("data", file3))) {
  download.file(paste0("https://raw.githubusercontent.com/lime-n/test/main/raster/", file3),
                file.path("data", file3), mode = "wb")
}

ebird_buff <- read_sf("data/shape.shp")
r <- raster("data/modis_mcd12q1_umd_2010.tif")

0
0

Please try and make reproducible examples. Reduce your pipe chains down to a minimum that creates the problem. Supply us with sample data where possible. Remove all superfluous code. That makes it so much easier for us.

Here is a simple data frame with two points taken from your Q:

> t_water = data.frame(X=c(-7158197,-7158197),Y=c(5814408,5812328))

Eventually your code does this - it creates an sf data frame and sets it to be coordinates in the epsg:4326 system:

> st_as_sf(t_water, coords = c("X", "Y"), crs = 4326)
Simple feature collection with 2 features and 0 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -7158197 ymin: 5812328 xmax: -7158197 ymax: 5814408
geographic CRS: WGS 84
                  geometry
1 POINT (-7158197 5814408)
2 POINT (-7158197 5812328)
> 

BUT epsg:4326 is degrees lat-long, and those numbers are clearly not degrees. If they are the same as the raster at that point then you don't need that step, you can do:

st_as_sf(t_water, coords = c("X", "Y"), crs = projection(landcover))

and then you also don't need the st_transform step since the X and Y are in that coordinate system, and you told st_as_sf that fact in that line.

It is so much easier to debug these things if you dont make long pipe chains. Save intermediate forms and then check them for correctness and consistency. Don't let garbage created in the middle of your pipe create garbage at the end.

4
  • You're right, I have tried what you've mentioned before, I think I might have mentioned this. Although, I will cut it down to provide detail of my problem. Given all of this, I'm still getting NA's for my raster. As for a reproducible example, it's fairly difficult without uploading an example raster (to download), and my dataset is fairly large.
    – Lime
    Mar 22 at 17:14
  • If you cant upload an example raster make some code that creates one, maybe with the same extent but fewer pixels, and fill it with 1:N. Should take two lines. Sometimes creating such data reveals the underlying misplaced assumption that is causing your problem in the first place.
    – Spacedman
    Mar 22 at 17:21
  • I have uploaded the code to download the files, I find recreating the data where this problem occurs bit tricky for my skills currently. Hopefully this is helpful
    – Lime
    Mar 22 at 17:39
  • I managed to figure out the problem earlier, though it was currently processing when I had this question up, such that it could be solved here instead of my 'back-up' plan. Essentially, I reprojected the habitat covariates that I extracted water values from to the same raster as r, then following the same code as above, and it provided me with the values I needed.
    – Lime
    Mar 22 at 23:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.