2

I am new-ish to using the stars package in R and I am having trouble with figuring out how to create a stars object C that has the values from object A but the extent of object B. Specifically, I have a map of average spring temperatures in Europe (Object A) and I want to crop it using a separate stars object containing presence of deciduous broadleaf forests (Object B).

Object A: https://i.stack.imgur.com/DQsZn.jpg

> CRU.SpringT.2009.2018_EU
stars object with 2 dimensions and 1 attribute
attribute(s):
                                   Min.  1st Qu.   Median     Mean  3rd Qu.     Max.   NA's
CRU.SpringT.2009.2018_EU.tif  -15.81895 1.201286 5.480992 4.979221 8.204463 17.55605 310479
dimension(s):
  from   to offset      delta refsys point values x/y
x    1 1440    -10  0.0416667 WGS 84 FALSE   NULL [x]
y    1  672     65 -0.0416667 WGS 84 FALSE   NULL [y]

Object B: https://i.stack.imgur.com/6dons.jpg

> Dec.BL_EU3
stars_proxy object with 1 attribute in 1 file(s):
$Consensus_reduced_class_3.tif
[1] "[...]/Consensus_reduced_class_3.tif"

dimension(s):
   from    to offset       delta                       refsys point values x/y
x 20401 27600   -180  0.00833333 +proj=longlat +datum=WGS8... FALSE   NULL [x]
y  3001  6360     90 -0.00833333 +proj=longlat +datum=WGS8... FALSE   NULL [y]
call_list:
[[1]]
x[i = i, drop = drop, crop = crop]
attr(,".Environment")
<environment: 0x000002a17a1ecac0>

[[2]]
x[i = i, drop = drop, crop = crop]
attr(,".Environment")
<environment: 0x000002a178d61dc0>

[[3]]
e1/e2
attr(,".Environment")
<environment: 0x000002a177987b10>

Both objects are cropped using the same bbox. The goal is for the resulting object (Object C) to have the extent of B but with the temperature values from A.

Dropbox links to tif files:

Object A: https://www.dropbox.com/s/lwvdxnis7k38e18/CRU.SpringT.2009.2018_EU.tif?dl=0

Object B: https://www.dropbox.com/s/uybxk40z853mu7a/EU%20Dec%20Broadleaf.tif?dl=0

Ben Lee
  • 35
  • 5
  • Hi @Ben Lee. Please provide the output summary of `Object A` and `Object B`. Without this information, it is difficult to help – lovalery Oct 20 '21 at 15:20
  • Hi @lovalery, I have edited my submission to include the output summaries. – Ben Lee Oct 20 '21 at 15:29
  • Not sure I understand. You say that `Both objects are cropped using the same bbox.` So A and B have the same extent. If you want to have the temperatures of A with the extent of B in an object C, that's like wanting to have the temperatures of A in the extent of A. And in this case it is enough to duplicate object_A like this `object_C <- object_A` Am I missing something? – lovalery Oct 20 '21 at 15:46
  • @lovalery Sorry, I am probably describing what I want using the wrong language. I want a stars object that looks like B (that is, it only has values where deciduous broadleaf forests exist), but, instead of 1, the values displayed are the temperature values from A. Right now B is just presence/absence of forests, but I want a continuous temperature value wherever presence = 1. – Ben Lee Oct 20 '21 at 16:06
  • O.K. thanks. It is much clearer for me now :-). Not easy without having the data at hand, but I will try to give you an answer to try (at least!) to put you on the right track :-) – lovalery Oct 20 '21 at 16:29

1 Answers1

1

That's it! It's much easier with the files :-) Please find below one possible solution to your problem. It works on my computer. I hope it will be the same for you :-)

I have deliberately detailed the code to make it easier for you to understand it. Of course, you can modify the form to make it more compact if needed.

PSEUDO Reprex !!

library(sf)
library(stars)
library(ggplot2)

A <- read_stars("CRU.SpringT.2009.2018_EU.tif")
B <- read_stars("consensus_full_class_3.tif")

# First step: Transform B to have the same extent and resolution of A
B_cropped_resample <- st_warp(B, A, use_gdal = TRUE)


# Second step: convert 0 to NA values in your object B (i.e. deciduous forests)
B_cropped_resample[B_cropped_resample == 0] <- NA # B should contain only NA and 1 values 


# Third step: convert B into `sf` object with 'points' geometry
B_cropped_resample_sf <- st_as_sf(B_cropped_resample, as_points = TRUE, na.rm = TRUE)


# Fourth step : extract values of object A with the object `B_cropped_resample_sf`
C <- st_extract(A, B_cropped_resample_sf) # C is a sf object which contains the 
                                          # values of A at the corresponding points

# Sixth step: convert 'sf' object back into stars 'object' and split 'sf' point 
# geometry into two 'stars' dimensions
C <- st_as_stars(C, name = attr(C, "CRU.SpringT.2009.2018_EU.tif", "geometry"))
C <- st_sfc2xy(C)


# Seventh step: plot the resulting 'stars' object which shows the temperatures (i.e.  
# values from object A) for each deciduous forest (i.e. locations from object B)
ggplot2::ggplot() + geom_stars(data = C) +
  coord_equal()

Please find the image obtained after having run the code:

enter image description here

lovalery
  • 4,524
  • 3
  • 14
  • 28
  • Hope this will help you. Please keep me informed of the result. Cheers. – lovalery Oct 20 '21 at 16:50
  • Hi @lovalery, this seems to be getting very close to what I want! However, the result is currently a list of points instead of a stars object (at least I think): `> C Simple feature collection with 6773164 features and 1 field Geometry type: POINT Dimension: XY Bounding box: xmin: -9.995833 ymin: 37.00417 xmax: 49.99583 ymax: 64.99583 CRS: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0` Is there a way to transform it back into a stars object? – Ben Lee Oct 20 '21 at 18:02
  • Yes of course :-) I will complete my answer. – lovalery Oct 20 '21 at 18:09
  • Edit made. Please keep me informed of the result! If everything is fine, please consider mark this answer as accepted. If not, feel free to tell me what is still going wrong. Cheers. – lovalery Oct 20 '21 at 18:14
  • Thanks, @lovalery, but now I have a new problem. I am hoping to plot this object using ggplot2 using `ggplot() + geom_stars(data = C)`, but when I try to do that with the new object I get an error `Error in check.length(gparname): 'gpar' element 'fontsize' must not be length 0` and a warning `In addition: Warning message: Removed 6773164 rows containing missing values (geom_sf)`. I think it might be because the object only has one dimension (see following comment) – Ben Lee Oct 20 '21 at 18:44
  • `> C stars object with 1 dimensions and 1 attribute attribute(s), summary of first 1e+05 cells: Min. 1st Qu. Median Mean 3rd Qu. Max. NA's CRU.SpringT.2009.2018_EU.tif -4.925655 -2.696862 -2.235813 -2.039258 -1.511543 3.714211 499 dimension(s): from to offset delta refsys point values geometry 1 6773164 NA NA WGS 84 TRUE POINT (11.8125 64.99583),...,POINT (49.99583 37.00417)` – Ben Lee Oct 20 '21 at 18:44
  • Could you provide me with the output summary of the sf object `B_sf`? – lovalery Oct 20 '21 at 18:50
  • `> B_sf Simple feature collection with 6773164 features and 1 field Geometry type: POINT` – Ben Lee Oct 20 '21 at 18:55
  • What is the name of the column (i.e. field) containing the values in the object `B_sf`? – lovalery Oct 20 '21 at 19:06
  • `Consensus_reduced_class_3.tif` – Ben Lee Oct 20 '21 at 19:27
  • No, you provide me with the name of your file. I just want the name of the column (i.e. field) containing the values in the object `B_sf`. You should have access to this column name when you display the output summary of the object as it shows the first 10 features. – lovalery Oct 20 '21 at 19:34
  • Hi @lovalery, I know that it is the name of the file, but for whatever reason it is also the name of the field. Please refer to this picture: https://i.stack.imgur.com/uREWx.jpg – Ben Lee Oct 20 '21 at 19:40
  • O.K. sorry. So, maybe try this when you convert the `sf` object to a `stars` object: `C <- st_as_stars(C, name = attr(C, "Consensus_reduced_class_3.tif"))` Keep me informed if it works better when you use `ggplot` – lovalery Oct 20 '21 at 19:47
  • Hi @lovalery, that solution did not work unfortunately and I am getting the same error message as before. – Ben Lee Oct 20 '21 at 23:53
  • Hello @Ben Lee, O.K. I made an edit in my answer above concerning the second step. The idea is to keep only the coordinates of the points in the `sf` object before extracting the values from the `stars` object. I keep my fingers crossed! Keep me informed of the result. – lovalery Oct 21 '21 at 09:40
  • Unfortunately this is still not working and I am getting the same error message as before when I try to plot it in ggplot. Maybe I should try using a different set of packages or something? Do you think this might be easier to do with raster instead? – Ben Lee Oct 21 '21 at 14:10
  • If possible, I think it would be more efficient if you provide me with your two already cropped objects/files (via a download platform like dropbox) so that I can test and come back to you with a proven solution. – lovalery Oct 21 '21 at 14:17
  • Hi @lovalery, thank you for your willingness to help. I have uploaded the two tifs to dropbox and edited my original question to link them at the end. Please let me know if you have difficulty downloading/accessing the files – Ben Lee Oct 21 '21 at 16:28
  • Thank you @Ben Lee. Well received. I will look into this and hope to get back to you soon with a solution :-) – lovalery Oct 21 '21 at 17:57
  • Hi @Ben Lee. I unfortunately come back with bad news. I tried to read the files (without spaces in their names) with `read_stars(x, proxy = FALSE)`. The file `CRU.SpringT.2009.2018_EU.tif` opens correctly but the `EU_Dec_Broadleaf.tif` file does not open (GDAL errors). I managed to read this last file with `terra::rast(EU_Dec_Broadleaf.tif)` but when converting to a `stars` object with the `st_as_stars()` function my pc crashed completely (Blue Screen of Death!!). – lovalery Oct 22 '21 at 11:57
  • I think this file is somewhat corrupted... On your side, do you manage to read this file by loading it into memory (i.e. `stars` object not a `stars proxy` object)? – lovalery Oct 22 '21 at 11:58
  • Hi @lovalery, my code for reading in that file looks like this: `Dec.BL <- stars::read_stars(paste0(fp.shp, "/Consensus_reduced_class_3.tif"))`. If the file I put in the dropbox is corrupted, maybe you will have better luck using the file straight from the source? I am using the third class from this website: https://www.earthenv.org/landcover – Ben Lee Oct 22 '21 at 12:26
  • O.K. thank you I will give it a try from the source. I will keep you informed of the result. – lovalery Oct 22 '21 at 12:34
  • Here is the solution (see above the answer completely rewritten). For the object B, I advise you to start from the source file because the file you sent me by the dropbox is weird (size much bigger than the file that can be loaded on the site you indicated...). In principle, everything should work correctly. I keep my fingers crossed :-). If everything works, please consider marking this answer as accepted. Cheers. – lovalery Oct 22 '21 at 16:36
  • This worked and is exactly what I wanted!!! Thank you @lovalery! – Ben Lee Oct 24 '21 at 14:57
  • Glad I could help you. Thank you very much for the validation and your vote. I wish you all the best with your work. Cheers. – lovalery Oct 24 '21 at 17:01