ropensci / MODIStsp

An "R" package for automatic download and preprocessing of MODIS Land Products Time Series
https://docs.ropensci.org/MODIStsp
GNU General Public License v3.0
155 stars 50 forks source link

Error in .local(.Object, ...) when trying to use MODIStsp_extract #236

Closed jschon637 closed 3 years ago

jschon637 commented 3 years ago

Issue description

I am unable to extract a time series from my Raster Stack, using MODIStsp_extract. I am struggling to figure out what I am doing wrong. I am sharing details about the objects I am using and my code. Any help to fix this will be greatly appreciated.

Also, my error is similar I believe to what appeared here: https://github.com/ropensci/MODIStsp/issues/173

Reproducible example


**Expected and actual behavior**
<!-- Provide here the full output of the provided example and describe what is going wrong- -->
setwd("C:/Users/jscho/OneDrive/NDVI")

#RasterStack file
infile    <- 'Africa_2016_200301-200501/VI_Monthly_1Km_v6/Time_Series/RData/Mixed/NDVI/MOD13A3_MYD13A3_NDVI_1_2003_1_2005_RData.RData' 

church<- st_read("Study_sites/study_sites.shp")
church<- st_transform(church,crs=4326)

######## #Set the start/end dates for extraction
startdate <- as.Date("2003-01-01")  
enddate   <- as.Date("2005-01-01")    

######## Load the RasterStack
inrts     <- get(load(infile)) 

######## Compute average
dataavg <- MODIStsp_extract(inrts,"Study_sites/study_sites.shp",
                            startdate, enddate, FUN = 'mean', na.rm = T)

### The output from the MODIStsp_extract step above is where I get the error, but below is more information
### about the church and inrts objects

> church
Simple feature collection with 47 features and 4 fields
geometry type:  POINT
dimension:      XYZ
bbox:           xmin: 37.46385 ymin: 11.52959 xmax: 38.07667 ymax: 11.98275
z_range:        zmin: 0 zmax: 0
geographic CRS: WGS 84
First 10 features:
   GISID_14 Sample      X      Y                      geometry
1      1025     25 37.770 37.770 POINT Z (37.76999 11.97073 0)
2      1028     59 37.751 37.751 POINT Z (37.75082 11.93424 0)
3      1032     24 37.741 37.741 POINT Z (37.74107 11.91171 0)
4      1034     60 37.702 37.702 POINT Z (37.70168 11.88704 0)
5      1037     61 37.672 37.672  POINT Z (37.67161 11.8413 0)
6      1054     62 37.997 37.997  POINT Z (37.99651 11.9146 0)
7      1074     63 38.007 38.007 POINT Z (38.00673 11.89767 0)
8      1075     64 38.001 38.001 POINT Z (38.00062 11.90513 0)
9      1076     37 38.006 38.006 POINT Z (38.00609 11.91433 0)
10     1077     39 38.028 38.028 POINT Z (38.02832 11.90295 0)

> inrts
class      : RasterStack 
dimensions : 8221, 8664, 71226744, 50  (nrow, ncol, ncell, nlayers)
resolution : 0.01025581, 0.01025546  (x, y)
extent     : -25.36055, 63.49575, -46.96973, 37.34041  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      : MOD13A3_NDVI_2003_001, MYD13A3_NDVI_2003_001, MOD13A3_NDVI_2003_032, MYD13A3_NDVI_2003_032, MOD13A3_NDVI_2003_060, MYD13A3_NDVI_2003_060, MOD13A3_NDVI_2003_091, MYD13A3_NDVI_2003_091, MOD13A3_NDVI_2003_121, MYD13A3_NDVI_2003_121, MOD13A3_NDVI_2003_152, MYD13A3_NDVI_2003_152, MOD13A3_NDVI_2003_182, MYD13A3_NDVI_2003_182, MOD13A3_NDVI_2003_213, ... 
time        : 2003-01-01 - 2005-01-01 (range)

I get the following output. I am specifically concerned about the line "Error in .local(.Object, ...) : "

although coordinates are longitude/latitude, st_intersection assumes that they are planar Error in .local(.Object, ...) : In addition: Warning message: attribute variables are assumed to be spatially constant throughout all geometries


**System information**
<!-- Provide here the output of the following R commands:
sessionInfo()
packageVersion("MODIStsp")
 -->

sessionInfo() R version 4.0.2 (2020-06-22) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale: [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages: [1] stats graphics grDevices utils datasets methods base

other attached packages: [1] cshapes_0.6 plyr_1.8.6 maptools_1.0-1 sf_0.9-3
[5] xts_0.12-0 zoo_1.8-8 raster_3.4-5 sp_1.4-1
[9] MODIStsp_2.0.5.9002

loaded via a namespace (and not attached): [1] Rcpp_1.0.5 xml2_1.3.2 magrittr_2.0.1 units_0.6-6
[5] lattice_0.20-41 R6_2.4.1 stringr_1.4.0 httr_1.4.2
[9] tools_4.0.2 rgdal_1.4-8 parallel_4.0.2 grid_4.0.2
[13] data.table_1.12.8 KernSmooth_2.23-17 e1071_1.7-3 DBI_1.1.0
[17] rgeos_0.5-3 remotes_2.3.0 class_7.3-17 yaml_2.2.1
[21] assertthat_0.2.1 gdalUtilities_1.1.2 bitops_1.0-6 codetools_0.2-16
[25] stringi_1.4.6 compiler_4.0.2 classInt_0.4-3 jsonlite_1.7.2
[29] foreign_0.8-80

packageVersion("MODIStsp") [1] ‘2.0.5.9002’

jschon637 commented 3 years ago

I'm not sure if this helps with my question, but the following code produces a similar error message:

raster::plot(inrts) Error in .local(.Object, ...) : Warning message: In graphics::par(old.par) : calling par(new=TRUE) with no plot

I also checked whether I could look at a subset of the RasterStack, and this is possible:

inrts[[1]] class : RasterLayer dimensions : 8221, 8664, 71226744 (nrow, ncol, ncell) resolution : 0.01025581, 0.01025546 (x, y) extent : -25.36055, 63.49575, -46.96973, 37.34041 (xmin, xmax, ymin, ymax) crs : +proj=longlat +datum=WGS84 +no_defs source : /home/jss5yf/Africa_2016/VI_Monthly_1Km_v6/NDVI/MOD13A3_NDVI_2003_001.tif names : MOD13A3_NDVI_2003_001 time : 2003-01-01

The source file name here comes from doing his download on my university's HPC.

jschon637 commented 3 years ago

Apologies, I stumbled on the solution. I had used MODIStsp to download the data on my university's HPC, and then downloaded the files onto my computer. Apparently, this messes up the RasterStack objects because they store summary information about the rasters, but not the raster data itself. So when I stack individual rasters myself on my computer, and then extract from that with MODIStsp_extract, everything works great. It looks like my problem came from misunderstanding RasterStack objects. I will close the issue now.