search

r-spatial / stars

Spatiotemporal Arrays, Raster and Vector Data Cubes
https://r-spatial.github.io/stars/
Apache License 2.0
563 stars 94 forks source link

Possible mistake in the use of sentinel bands to calculate NDVI. Vignettes stars::stars2 (2. stars proxy objects) #336

Closed alexyshr closed 4 years ago

alexyshr commented 4 years ago

According to this, the description for the 10 meters bands of sentinel 2 are:

  1. B2 - Blue
  2. B3 - Green
  3. B4 - Red
  4. B8 - Visible and Near Infrared (VNIR)

So the function to calculate NDVI ((NIR – RED) / (NIR + RED)) in Vignette Stars Proxy Objects need to be changed from: ndvi = function(x) (x[4] - x[1])/(x[4] + x[1]) to: ndvi = function(x) (x[4] - x[3])/(x[4] + x[3])

I am using gdalUtils::gdalinfo to check the content of the Sentinel 2 file present in starsdatapackage (not available in R version 4.0.2)

library(starsdata)
library(stars)
#> Warning: package 'stars' was built under R version 3.6.3
#> Loading required package: abind
#> Loading required package: sf
#> Warning: package 'sf' was built under R version 3.6.3
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(gdalUtils)
#> Warning: package 'gdalUtils' was built under R version 3.6.3
#> 
#> Attaching package: 'gdalUtils'
#> The following object is masked from 'package:sf':
#> 
#>     gdal_rasterize
granule <- system.file("sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip", package = "starsdata")
s2 <- paste0("SENTINEL2_L1C:/vsizip/", granule, "/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/MTD_MSIL1C.xml:10m:EPSG_32632")
(p <- read_stars(s2, proxy = TRUE))
#> stars_proxy object with 1 attribute in file:
#> $`MTD_MSIL1C.xml:10m:EPSG_32632`
#> [1] "[...]/MTD_MSIL1C.xml:10m:EPSG_32632"
#> 
#> dimension(s):
#>      from    to offset delta                refsys point values    
#> x       1 10980  3e+05    10 WGS 84 / UTM zone 32N    NA   NULL [x]
#> y       1 10980  6e+06   -10 WGS 84 / UTM zone 32N    NA   NULL [y]
#> band    1     4     NA    NA                    NA    NA   NULL

gdalinfo(granule)
#>  [1] "Driver: SENTINEL2/Sentinel 2"                                                                                                                                                                                                                                                                                                                                                                                 
#>  [2] "Files: E:/Documents/R/win-library/3.6/starsdata/sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip"                                                                                                                                                                                                                                                                                    
#>  [3] "Size is 512, 512"                                                                                                                                                                                                                                                                                                                                                                                             
#>  [4] "Coordinate System is `'"                                                                                                                                                                                                                                                                                                                                                                                      
#>  [5] "Metadata:"                                                                                                                                                                                                                                                                                                                                                                                                    
#>  [6] "  CLOUD_COVERAGE_ASSESSMENT=66.2526"                                                                                                                                                                                                                                                                                                                                                                          
#>  [7] "  DATATAKE_1_DATATAKE_SENSING_START=2018-02-20T10:50:51.026Z"                                                                                                                                                                                                                                                                                                                                                 
#>  [8] "  DATATAKE_1_DATATAKE_TYPE=INS-NOBS"                                                                                                                                                                                                                                                                                                                                                                          
#>  [9] "  DATATAKE_1_ID=GS2A_20180220T105051_013919_N02.06"                                                                                                                                                                                                                                                                                                                                                           
#> [10] "  DATATAKE_1_SENSING_ORBIT_DIRECTION=DESCENDING"                                                                                                                                                                                                                                                                                                                                                              
#> [11] "  DATATAKE_1_SENSING_ORBIT_NUMBER=51"                                                                                                                                                                                                                                                                                                                                                                         
#> [12] "  DATATAKE_1_SPACECRAFT_NAME=Sentinel-2A"                                                                                                                                                                                                                                                                                                                                                                     
#> [13] "  DEGRADED_ANC_DATA_PERCENTAGE=0.0"                                                                                                                                                                                                                                                                                                                                                                           
#> [14] "  DEGRADED_MSI_DATA_PERCENTAGE=0"                                                                                                                                                                                                                                                                                                                                                                             
#> [15] "  FOOTPRINT=POLYGON((7.631851463227346 53.750197934337095, 7.597710703358927 53.69050276665881, 7.516336043243608 53.547183540361466, 7.435866385199494 53.403802271794845, 7.355637002948162 53.26039397065163, 7.29250760569914 53.14699997423984, 6.010910777116657 53.123645037009204, 5.940406316632895 54.109206448119856, 7.619256720590609 54.140186775460364, 7.631851463227346 53.750197934337095))"
#> [16] "  FORMAT_CORRECTNESS=PASSED"                                                                                                                                                                                                                                                                                                                                                                                  
#> [17] "  GENERAL_QUALITY=PASSED"                                                                                                                                                                                                                                                                                                                                                                                     
#> [18] "  GENERATION_TIME=2018-02-21T13:40:37.000000Z"                                                                                                                                                                                                                                                                                                                                                                
#> [19] "  GEOMETRIC_QUALITY=PASSED"                                                                                                                                                                                                                                                                                                                                                                                   
#> [20] "  PREVIEW_GEO_INFO=Not applicable"                                                                                                                                                                                                                                                                                                                                                                            
#> [21] "  PREVIEW_IMAGE_URL=Not applicable"                                                                                                                                                                                                                                                                                                                                                                           
#> [22] "  PROCESSING_BASELINE=02.06"                                                                                                                                                                                                                                                                                                                                                                                  
#> [23] "  PROCESSING_LEVEL=Level-1C"                                                                                                                                                                                                                                                                                                                                                                                  
#> [24] "  PRODUCT_START_TIME=2018-02-20T10:50:51.026Z"                                                                                                                                                                                                                                                                                                                                                                
#> [25] "  PRODUCT_STOP_TIME=2018-02-20T10:50:51.026Z"                                                                                                                                                                                                                                                                                                                                                                 
#> [26] "  PRODUCT_TYPE=S2MSI1C"                                                                                                                                                                                                                                                                                                                                                                                       
#> [27] "  PRODUCT_URI=S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE"                                                                                                                                                                                                                                                                                                                              
#> [28] "  QUANTIFICATION_VALUE=10000"                                                                                                                                                                                                                                                                                                                                                                                 
#> [29] "  RADIOMETRIC_QUALITY=PASSED"                                                                                                                                                                                                                                                                                                                                                                                 
#> [30] "  REFLECTANCE_CONVERSION_U=1.02459181087746"                                                                                                                                                                                                                                                                                                                                                                  
#> [31] "  SENSOR_QUALITY=PASSED"                                                                                                                                                                                                                                                                                                                                                                                      
#> [32] "  SPECIAL_VALUE_NODATA=0"                                                                                                                                                                                                                                                                                                                                                                                     
#> [33] "  SPECIAL_VALUE_SATURATED=65535"                                                                                                                                                                                                                                                                                                                                                                              
#> [34] "Subdatasets:"                                                                                                                                                                                                                                                                                                                                                                                                 
#> [35] "  SUBDATASET_1_NAME=SENTINEL2_L1C:/vsizip/E:/Documents/R/win-library/3.6/starsdata/sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/MTD_MSIL1C.xml:10m:EPSG_32632"                                                                                                                                                 
#> [36] "  SUBDATASET_1_DESC=Bands B2, B3, B4, B8 with 10m resolution, UTM 32N"                                                                                                                                                                                                                                                                                                                                        
#> [37] "  SUBDATASET_2_NAME=SENTINEL2_L1C:/vsizip/E:/Documents/R/win-library/3.6/starsdata/sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/MTD_MSIL1C.xml:20m:EPSG_32632"                                                                                                                                                 
#> [38] "  SUBDATASET_2_DESC=Bands B5, B6, B7, B8A, B11, B12 with 20m resolution, UTM 32N"                                                                                                                                                                                                                                                                                                                             
#> [39] "  SUBDATASET_3_NAME=SENTINEL2_L1C:/vsizip/E:/Documents/R/win-library/3.6/starsdata/sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/MTD_MSIL1C.xml:60m:EPSG_32632"                                                                                                                                                 
#> [40] "  SUBDATASET_3_DESC=Bands B1, B9, B10 with 60m resolution, UTM 32N"                                                                                                                                                                                                                                                                                                                                           
#> [41] "  SUBDATASET_4_NAME=SENTINEL2_L1C:/vsizip/E:/Documents/R/win-library/3.6/starsdata/sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/MTD_MSIL1C.xml:TCI:EPSG_32632"                                                                                                                                                 
#> [42] "  SUBDATASET_4_DESC=True color image, UTM 32N"                                                                                                                                                                                                                                                                                                                                                                
#> [43] "Corner Coordinates:"                                                                                                                                                                                                                                                                                                                                                                                          
#> [44] "Upper Left  (    0.0,    0.0)"                                                                                                                                                                                                                                                                                                                                                                                
#> [45] "Lower Left  (    0.0,  512.0)"                                                                                                                                                                                                                                                                                                                                                                                
#> [46] "Upper Right (  512.0,    0.0)"                                                                                                                                                                                                                                                                                                                                                                                
#> [47] "Lower Right (  512.0,  512.0)"                                                                                                                                                                                                                                                                                                                                                                                
#> [48] "Center      (  256.0,  256.0)"

# According  with gdalinfo (MTD_MSIL1C.xml:10m:EPSG_32632):
# [35] "  SUBDATASET_1_NAME=SENTINEL2_L1C:/vsizip/E:/Documents/R/win-library/3.6/starsdata/sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/MTD_MSIL1C.xml:10m:EPSG_32632"
# [36] "  SUBDATASET_1_DESC=Bands B2, B3, B4, B8 with 10m resolution, UTM 32N"

# The sub-dataset one is composed of bands: B2, B3, B4, B8, so the four bands in stars 'p' are:
# 1: B2
# 2: B3
# 3: B4
# 4: B8

# And according to <https://gisgeography.com/sentinel-2-bands-combinations/#:~:text=Sentinel%2D2%20carries%20the%20Multispectral,to%2060%2Dmeter%20pixel%20size.&text=Next%2C%20its%20red%20edge%20(B5,sampling%20distance%20of%2020%20meters.>
# 1: B2: Blue
# 2: B3: Green
# 3: B4: Red
# 4: B8: Visible and Near Infrared (VNIR)

# So next part in the Vignette stars::stars2 (2. stars proxy objects):
# "S2 10m: band 4: near infrared, band 1: red."
# need to be changed to:
# "S2 10m: band 4: near infrared, band 3: red"
# The band 1 is blue not red. The red band is 3
# And the ndvi function need to be changed to:
# ndvi = function(x) (x[4] - x[3])/(x[4] + x[3])
sessionInfo()
#> R version 3.6.2 (2019-12-12)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 18363)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.1252 
#> [2] LC_CTYPE=English_United States.1252   
#> [3] LC_MONETARY=English_United States.1252
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.1252    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> loaded via a namespace (and not attached):
#>  [1] compiler_3.6.2  magrittr_1.5    tools_3.6.2     htmltools_0.4.0
#>  [5] yaml_2.2.1      Rcpp_1.0.5      stringi_1.4.6   rmarkdown_2.2  
#>  [9] highr_0.8       knitr_1.28      stringr_1.4.0   xfun_0.14      
#> [13] digest_0.6.25   rlang_0.4.7     evaluate_0.14
edzer commented 4 years ago

With stars installed from github, I see

> st_dimensions(p)$band$values
[1] "B4" "B3" "B2" "B8"

Also:

> gdal_utils("info", p)
Error in CPL_gdalinfo(source, options, oo) : 
  Not compatible with STRSXP: [type=list].
> gdal_utils("info", s2)
Driver: SENTINEL2/Sentinel 2
Files: /vsizip//home/edzer/R/x86_64-pc-linux-gnu-library/4.0/starsdata/sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/MTD_MSIL1C.xml
       /vsizip//home/edzer/R/x86_64-pc-linux-gnu-library/4.0/starsdata/sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/GRANULE/L1C_T32ULE_A013919_20180220T105539/MTD_TL.xml
       /vsizip//home/edzer/R/x86_64-pc-linux-gnu-library/4.0/starsdata/sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/GRANULE/L1C_T32ULE_A013919_20180220T105539/IMG_DATA/T32ULE_20180220T105051_B04.jp2
       /vsizip//home/edzer/R/x86_64-pc-linux-gnu-library/4.0/starsdata/sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/GRANULE/L1C_T32ULE_A013919_20180220T105539/IMG_DATA/T32ULE_20180220T105051_B03.jp2
       /vsizip//home/edzer/R/x86_64-pc-linux-gnu-library/4.0/starsdata/sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/GRANULE/L1C_T32ULE_A013919_20180220T105539/IMG_DATA/T32ULE_20180220T105051_B02.jp2
       /vsizip//home/edzer/R/x86_64-pc-linux-gnu-library/4.0/starsdata/sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/GRANULE/L1C_T32ULE_A013919_20180220T105539/IMG_DATA/T32ULE_20180220T105051_B08.jp2
Size is 10980, 10980
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 32N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 32N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    ID["EPSG",32632]]
Data axis to CRS axis mapping: 1,2
Origin = (300000.000000000000000,6000000.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  CLOUD_COVERAGE_ASSESSMENT=66.2526
  DATATAKE_1_DATATAKE_SENSING_START=2018-02-20T10:50:51.026Z
  DATATAKE_1_DATATAKE_TYPE=INS-NOBS
  DATATAKE_1_ID=GS2A_20180220T105051_013919_N02.06
  DATATAKE_1_SENSING_ORBIT_DIRECTION=DESCENDING
  DATATAKE_1_SENSING_ORBIT_NUMBER=51
  DATATAKE_1_SPACECRAFT_NAME=Sentinel-2A
  DEGRADED_ANC_DATA_PERCENTAGE=0.0
  DEGRADED_MSI_DATA_PERCENTAGE=0
  FORMAT_CORRECTNESS=PASSED
  GENERAL_QUALITY=PASSED
  GENERATION_TIME=2018-02-21T13:40:37.000000Z
  GEOMETRIC_QUALITY=PASSED
  PREVIEW_GEO_INFO=Not applicable
  PREVIEW_IMAGE_URL=Not applicable
  PROCESSING_BASELINE=02.06
  PROCESSING_LEVEL=Level-1C
  PRODUCT_START_TIME=2018-02-20T10:50:51.026Z
  PRODUCT_STOP_TIME=2018-02-20T10:50:51.026Z
  PRODUCT_TYPE=S2MSI1C
  PRODUCT_URI=S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE
  QUANTIFICATION_VALUE=10000
  RADIOMETRIC_QUALITY=PASSED
  REFLECTANCE_CONVERSION_U=1.02459181087746
  SENSOR_QUALITY=PASSED
  SPECIAL_VALUE_NODATA=0
  SPECIAL_VALUE_SATURATED=65535
Image Structure Metadata:
  COMPRESSION=JPEG2000
Corner Coordinates:
Upper Left  (  300000.000, 6000000.000) (  5d56'25.46"E, 54d 6'33.14"N)
Lower Left  (  300000.000, 5890200.000) (  6d 0'39.28"E, 53d 7'25.12"N)
Upper Right (  409800.000, 6000000.000) (  7d37' 9.32"E, 54d 8'24.67"N)
Lower Right (  409800.000, 5890200.000) (  7d39' 4.03"E, 53d 9'12.73"N)
Center      (  354900.000, 5945100.000) (  6d48'19.52"E, 53d38' 4.35"N)
Band 1 Block=128x128 Type=UInt16, ColorInterp=Red
  Description = B4, central wavelength 665 nm
  Overviews: 5490x5490, 2745x2745, 1373x1373, 687x687
  Metadata:
    BANDNAME=B4
    BANDWIDTH=30
    BANDWIDTH_UNIT=nm
    SOLAR_IRRADIANCE=1512.06
    SOLAR_IRRADIANCE_UNIT=W/m2/um
    WAVELENGTH=665
    WAVELENGTH_UNIT=nm
  Image Structure Metadata:
    NBITS=15
Band 2 Block=128x128 Type=UInt16, ColorInterp=Green
  Description = B3, central wavelength 560 nm
  Overviews: 5490x5490, 2745x2745, 1373x1373, 687x687
  Metadata:
    BANDNAME=B3
    BANDWIDTH=35
    BANDWIDTH_UNIT=nm
    SOLAR_IRRADIANCE=1823.24
    SOLAR_IRRADIANCE_UNIT=W/m2/um
    WAVELENGTH=560
    WAVELENGTH_UNIT=nm
  Image Structure Metadata:
    NBITS=15
Band 3 Block=128x128 Type=UInt16, ColorInterp=Blue
  Description = B2, central wavelength 490 nm
  Overviews: 5490x5490, 2745x2745, 1373x1373, 687x687
  Metadata:
    BANDNAME=B2
    BANDWIDTH=65
    BANDWIDTH_UNIT=nm
    SOLAR_IRRADIANCE=1959.72
    SOLAR_IRRADIANCE_UNIT=W/m2/um
    WAVELENGTH=490
    WAVELENGTH_UNIT=nm
  Image Structure Metadata:
    NBITS=15
Band 4 Block=128x128 Type=UInt16, ColorInterp=Undefined
  Description = B8, central wavelength 842 nm
  Overviews: 5490x5490, 2745x2745, 1373x1373, 687x687
  Metadata:
    BANDNAME=B8
    BANDWIDTH=115
    BANDWIDTH_UNIT=nm
    SOLAR_IRRADIANCE=1041.63
    SOLAR_IRRADIANCE_UNIT=W/m2/um
    WAVELENGTH=842
    WAVELENGTH_UNIT=nm
  Image Structure Metadata:
    NBITS=15

It would be nice if we could address the bands by name in the ndvi function, rather than by number.