Harmonized Landsat Sentinel HLS CRS problem

[s6pubrue]

Hello, I want to use HLS data for my master thesis and I'm currently following this tutorial: https://github.com/nasa/HLS-Data-Resources/blob/main/r/HLS_Tutorial.Rmd
In doing so, I come across a problem with this part here:


{r, warning=FALSE, results= "hide"}
red_stack <- nir_stack <- fmask_stack <- date_list <- list()
# Add progress bar
pb = txtProgressBar(min = 0, max = length(search_df$band), initial = 0, style = 3)
l <- m <- n <- 0
for (row in seq(length(search_df$band))){
setTxtProgressBar(pb,row)
if (search_df$band[row] == 'B04'){
l = l+1
red <- terra::rast(paste0('/vsicurl/', search_df$Asset_Link[row]))
red_crop <- terra::mask(terra::crop(red, terra::ext(cropland_geojson_utm)), cropland_geojson_utm)
red_stack[[l]] <- red_crop

doy_time = strsplit(sources(red), "[.]")[[1]][14]
doy = substr(doy_time, 1, 7)
date <- as.Date(as.integer(substr(doy,5,7)), origin = paste0(substr(doy,1,4), "-01-01"))

if (strsplit(sources(red), "[.]")[[1]][12] == 'S30'){
date_list[[l]] <- paste0('S',as.character(date))
}else{
date_list[[l]] <- paste0('L',as.character(date))
}

rm (red, red_crop)
}else if (search_df$band[row] == 'Fmask'){
m = m+1
fmask <- terra::rast(paste0('/vsicurl/', search_df$Asset_Link[row]))
fmask_crop <- terra::mask (terra::crop(fmask, terra::ext(cropland_geojson_utm)), cropland_geojson_utm)
fmask_stack[[m]] <- fmask_crop
rm(fmask, fmask_crop)
}else{
n = n+1
nir <- terra::rast(paste0('/vsicurl/', search_df$Asset_Link[row]))
nir_crop <- terra::mask (terra::crop(nir, terra::ext(cropland_geojson_utm)), cropland_geojson_utm)
nir_stack[[n]] <- nir_crop
rm(nir, nir_crop)
}
}
close(pb)


When this has run through, I am shown:


There were 21 warnings (display with warnings())
> close(pb)

> warnings()

1: [mask] CRS do not match
2: [mask] CRS do not match
3: [mask] CRS do not match
4: [mask] CRS do not match
5: [mask] CRS do not match
6: [mask] CRS do not match
7: [mask] CRS do not match
8: [mask] CRS do not match
9: [mask] CRS do not match
10: [mask] CRS do not match
11: [mask] CRS do not match
12: [mask] CRS do not match
13: [mask] CRS do not match
14: [mask] CRS do not match
15: [mask] CRS do not match
16: [mask] CRS do not match
17: [mask] CRS do not match
18: [mask] CRS do not match
19: [mask] CRS do not match
20: [mask] CRS do not match
21: [mask] CRS do not match

When I then run nir_stack, it shows:

[[1]]
class : SpatRaster
dimensions : 117, 178, 1 (nrow, ncol, nlyr)
resolution : 30, 30 (x, y)
extent : 577200, 582540, 4416150, 4419660 (xmin, xmax, ymin, ymax)
coord. ref. : UTM Zone 10, Northern Hemisphere
source(s) : memory
varname : HLS.S30.T10TEK.2021214T184919.v2.0.B8A
name : NIR_Narrow
min value : 0.0371
max value : 0.4857

[[2]]
class : SpatRaster
dimensions : 117, 178, 1 (nrow, ncol, nlyr)
resolution : 30, 30 (x, y)
extent : 577200, 582540, 4416150, 4419660 (xmin, xmax, ymin, ymax)
coord. ref. : UTM Zone 10, Northern Hemisphere
source(s) : memory
varname : HLS.S30.T10TEK.2021217T185919.v2.0.B8A
name : NIR_Narrow
min value : 0.0317
max value : 0.4847

[[3]]
class : SpatRaster
dimensions : 117, 178, 1 (nrow, ncol, nlyr)
resolution : 30, 30 (x, y)
extent : 577200, 582540, 4416150, 4419660 (xmin, xmax, ymin, ymax)
coord. ref. : UTM Zone 10, Northern Hemisphere
source(s) : memory
varname : HLS.S30.T10TEK.2021219T184921.v2.0.B8A
name : NIR_Narrow
min value : 0.0992
max value : 0.4219

[[4]]
class : SpatRaster
dimensions : 117, 178, 1 (nrow, ncol, nlyr)
resolution : 30, 30 (x, y)
extent : 577200, 582540, 4416150, 4419660 (xmin, xmax, ymin, ymax)
coord. ref. : UTM Zone 10, Northern Hemisphere
source(s) : memory
name : NIR
min value : 1112
max value : 3984

In total, search_df has 72 entries after cloud filtering. So these 4 stacks definately do not include all found images. When I then run the following part, I get an index out of range error:

ndvi_stack <- list()
for (i in 1:length(nir_stack)){
# Calculate NDVI
ndvi_stack[[i]] <- calculate_NDVI(nir_stack[[i]], red_stack[[i]])
# Exclude the Inf and -Inf from NDVI
ndvi_stack[[i]][ndvi_stack[[i]] == Inf] <- NA
ndvi_stack[[i]][ndvi_stack[[i]] == -Inf] <- NA
names(ndvi_stack[[i]]) <- date_list[[i]]
}
ndvi_stacks <- terra::rast(ndvi_stack) # Create a stack of NDVI


I then tested around a bit by inserting different numbers in search_df$Asset_Link[1] and checking the CRS. When executing

> coordinate_reference <- terra::rast(paste0('/vsicurl/',search_df$Asset_Link[1])) # Extract the CRS from one of the assets
> crs(coordinate_reference)

the following is displayed:

[1] "PROJCRS[\"UTM Zone 10, Northern Hemisphere\",\n BASEGEOGCRS[\"Unknown datum based upon the WGS 84 ellipsoid\",\n DATUM[\"Not specified (based on WGS 84 spheroid)\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",7030]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]]],\n CONVERSION[\"Transverse Mercator\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-123,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"easting\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n AXIS[\"northing\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]]"

But if I insert a 68 (I just randomly tested some numbers) instead of the 1 and execute this

> coordinate_reference <- terra::rast(paste0('/vsicurl/',search_df$Asset_Link[68])) # Extract the CRS from one of the assets
> crs(coordinate_reference)

then this is displayed:

[1] "PROJCRS[\"WGS 84 / UTM zone 10N\",\n BASEGEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]],\n CONVERSION[\"UTM zone 10N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-123,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n AREA[\"Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK).\"],\n BBOX[0,-126,84,-120]],\n ID[\"EPSG\",32610]]"

So both coordinate systems are UTM zone 10 North and are also both based on WGS 84, but they are slightly different. Could this be the reason for the warning messages? Was there possibly a change in the CRS? I followed all the steps in the tutorial exactly and just ran the code chunks, so the error should be reproducible. I would really appreciate your help as I urgently need the data soon. Thank you very much!

[LP DAAC-EDL - dgolon]

Hi @s6pubrue We are talking with our developers on this and will report back when we have an answer.

[LP DAAC-EDL - kyrad]

Hello @s6pubrue, are you able to send us your area of interest, time, and product you are using? We would like to try to recreate the error you are getting. Thanks

[LP DAAC-EDL - dgolon]

Hello, after discussing with our developers, this is expected. The descriptions in the CRS are different between HLS S30 and HLS L30. The CRS is the same, the WKT string is different, but the CRS is what you should focus on. The tutorial will throw a warning, but you can ignore it. Thanks -- Danielle