library(stars)

## Loading required package: abind

## Loading required package: sf

## Linking to GEOS 3.7.0, GDAL 2.4.0, PROJ 5.2.0

# Read 1km MODIS grid with sinusoidal CRS
r = read_stars("/media/qnap/Data/NDVI/US_Seasonal_NDVI_Landsat_120m_winter_2000_2018_summer_2000_2018/template.tif")

# Subset
r = r[, 2000:2050, 2000:2050, ]
r

## stars object with 2 dimensions and 1 attribute
## attribute(s):
##  template.tif     
##  Min.   :7859440  
##  1st Qu.:7904650  
##  Median :7951743  
##  Mean   :7950731  
##  3rd Qu.:7997154  
##  Max.   :8038642  
## dimension(s):
##   from   to    offset    delta                       refsys point values
## x 2000 2050 -11123212  926.625 +proj=sinu +lon_0=0 +x_0=... FALSE   NULL
## y 2000 2050   5493036 -926.625 +proj=sinu +lon_0=0 +x_0=... FALSE   NULL
##      
## x [x]
## y [y]

plot(r, main = "Sinusoidal - regular grid")

# Reproject - 4326
r1 = st_transform(r, 4326)
r1

## stars object with 2 dimensions and 1 attribute
## attribute(s):
##  template.tif     
##  Min.   :7859440  
##  1st Qu.:7904650  
##  Median :7951743  
##  Mean   :7950731  
##  3rd Qu.:7997154  
##  Max.   :8038642  
## dimension(s):
##   from   to offset delta                       refsys point
## x 2000 2050     NA    NA +proj=longlat +datum=WGS8... FALSE
## y 2000 2050     NA    NA +proj=longlat +datum=WGS8... FALSE
##                          values    
## x [51x51] -99.1145,...,-98.1628 [x]
## y   [51x51] 32.3208,...,32.7375 [y]
## curvilinear grid

plot(r1, main = "4326 - Curvilinear grid")

# Reproject - 2163
r2 = st_transform(r, 2163)
r2

## stars object with 2 dimensions and 1 attribute
## attribute(s):
##  template.tif     
##  Min.   :7859440  
##  1st Qu.:7904650  
##  Median :7951743  
##  Mean   :7950731  
##  3rd Qu.:7997154  
##  Max.   :8038642  
## dimension(s):
##   from   to offset delta                       refsys point
## x 2000 2050     NA    NA +proj=laea +lat_0=45 +lon... FALSE
## y 2000 2050     NA    NA +proj=laea +lat_0=45 +lon... FALSE
##                          values    
## x    [51x51] 83300.2,...,173685 [x]
## y [51x51] -1405988,...,-1359879 [y]
## curvilinear grid

plot(r2, main = "2163 - Curvilinear grid")