Comparison of WIP MARS pipeline to ECMWF-AWS NRT pipeline

Processed Ensemble Mean COGs

Intro

This document is a comparison of the MARS processed COGS and the AWS processed COGS.

  • ECMWF has been pushing 0.4 degree NRT seasonal forecast gribs to our AWS bucket on the 5th of each month (starting 2024). These gribs are processed into COGS by taking the ensemble mean.
  • We want to compare the above w/ seasonal forecast ensemble mean COGS from MARS catalogue. These COGS have been previously produced in an adhoc manner for AA projects such as Dry Corridor, Afghanistan, Ethiopia. But here we want to compare with the new WIP Global MARS Seasonal Forecast Ensemble Mean COG pipeline

Just recently I processed the global AWS grib file into ensemble mean cogs for July 2024 and pushed to our Azure Blob. For convenience, I will compare that to the MARS ensemble mean COGS.

Code 
box::use(blob = ../R/utils)
box::use(AzureStor)
box::use(dplyr[...])
box::use(lubridate[...])
box::use(stringr[...])
box::use(sf[...])
box::use(terra[...])
box::use(rnaturalearth[...])

Sys.setenv(AZURE_SAS = Sys.getenv("DSCI_AZ_SAS_DEV"))
Sys.setenv(AZURE_STORAGE_ACCOUNT = Sys.getenv("DSCI_AZ_STORAGE_ACCOUNT"))

extract_date <- function(x){
    as_date(str_extract(x,"\\d{4}-\\d{2}-\\d{2}"))
}


# helper function for aws bucket
get_formatted_bucket_df <-  function(bucket=Sys.getenv("BUCKET_NAME")){
    aws.s3::get_bucket_df(bucket,max = 10000) |> 
        dplyr::as_tibble() |> 
        dplyr::mutate(
            date= lubridate::as_date(LastModified)
        ) |> 
        tidyr::separate(LastModified,into = c("date","time"), sep = "T") %>% 
        dplyr::mutate(
            time= stringr::str_remove(time,"\\.000Z"),
            time = stringr::str_replace_all(time,":","-"),
            filename = paste0(date,"_",time,".grib"),
            date= lubridate::as_date(date),
            size_mb = as.numeric(Size) * 1e-6
        )
}

Load blobs

  • At current time we know that the global COGS produced from MARS GRIBS are in the blob for the year 2000. Therefore we will compare that to the COGS produced from ensemble means of AWS GRIBS from July 2024 to those produced from MARS GRIBS July 2000.
  • Note not currently possible to do a direct comparison of the same date COGS as MARS catalogue has long publishing delay need to find citation, whereas we have only recieved the NRT data on our AWS bucket since 2024.
  • As the data was processed slightly different we will just wrangle the MARS derived bands a bit to match the others for the purpose of this doc.
Code 
pc <- blob$load_proj_contatiners()

az_prefix <- "/vsiaz/global/"
aws_cog_dir = "raster/cogs"
mars_cog_dir = "mars/processed"


l_contents_aws <- AzureStor$list_blobs(
    container = pc$GLOBAL_CONT,
    prefix = paste0(aws_cog_dir,"/202407-ECMWF")
)

l_contents_mars <-  AzureStor$list_blobs(
    container = pc$GLOBAL_CONT,
    dir = mars_cog_dir
)

l_contents_mars <- l_contents_mars |> 
    mutate(
        date= extract_date(name),
        virt_path = paste0(az_prefix,name)
    ) 

# let just a random mars pub date collection to compare to the AWS
mars_urls <- l_contents_mars |> 
    filter(
        # data from 2000 is the new global data
        date == "2000-07-01"
    ) |> 
    pull(virt_path)

aws_url <-  paste0(az_prefix, l_contents_aws$name)
r_aws <- rast(aws_url)
r_mars <- rast(mars_urls)


# adjust r_mars names so they match r_aws for easier comparison
bnames_modified_mars <- 
    str_remove_all(
        basename(sources(r_mars)),
        "^seas5_mars_tprate_em_i|lt|\\.tif$"
    ) |> 
    str_replace_all("_",".")

names(r_mars) <- bnames_modified_mars

Check Consistency/Alignment

Both r_mars and r_aws produce the same warning message.

Code 
plot(r_mars[[1]])
Warning: [is.lonlat] coordinates are out of range for lon/lat

Code 
plot(r_aws[[1]])
Warning: [is.lonlat] coordinates are out of range for lon/lat

If the grids & extents are aligned we should be able to merge/stack them. The call below does produce a warning: [rast] CRS do not match. But nonetheless still allows the stacking of bands.

Code 
rast(
    list(
        r_mars[[1]],
        r_aws[[1]]
    )
)
Warning: [rast] CRS do not match
class       : SpatRaster 
dimensions  : 451, 900, 2  (nrow, ncol, nlyr)
resolution  : 0.4, 0.4  (x, y)
extent      : -180.2, 179.8, -90.2, 90.2  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : seas5_mars_tprate_em_i2000-07-01_lt0.tif  
              202407-ECMWF_SEAS-V_i_a-ensemble_mean.tif  
names       : 2000-07-01.0, 2024-07-01.0 
min values  :           ? ,        0.000 
max values  :           ? ,     1741.675

CRS exploration

It looks like the r_aws CRS is “imported from GRIB file” and not independently set, whereas the MARS one has actually been set.

Code 
st_crs(r_aws[[1]])
Coordinate Reference System:
  User input: Coordinate System imported from GRIB file 
  wkt:
GEOGCRS["Coordinate System imported from GRIB file",
    DATUM["unnamed",
        ELLIPSOID["Sphere",6367470,0,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122]]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]
Code 
st_crs(r_mars[[1]])
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

Adjusting AWS derived COG CRS

Two ways we can adjust the CRS of the AWS COG so it matches MARS COG

  1. we can project - this would be used if the CRS was in fact different and we wanted to project to new CRS
  2. we can simply set the CRS - this could be used if the coordinates are in fact correct for WGS84, but for some reason the registered CRS just isn’t aligning
Code 
# reproject
r_aws_proj <- terra::project(r_aws[[1]], crs(r_mars[[1]]))
Code 
# set
r_aws_copy <- deepcopy(r_aws)
crs(r_aws_copy) <-  crs(r_mars)

When trying to merge bands of reprojected AWS derived w/ MARS-derived, we get: Error: [rast] extents do not match. This error combined w/ the previous warning indicates that the grids were actually aligned previously, and this re-projections incorrectly changed the coordinates of the r_aws

Code 
rast(
    list(
        r_mars[[1]],
        r_aws_proj
    )
)
Error: [rast] extents do not match

This is then further confirmed by the fact that the r_aws_copy does not produce the error

Code 
rast(
    list(
        r_mars[[1]],
        r_aws_copy[[1]]
    )
)
class       : SpatRaster 
dimensions  : 451, 900, 2  (nrow, ncol, nlyr)
resolution  : 0.4, 0.4  (x, y)
extent      : -180.2, 179.8, -90.2, 90.2  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : seas5_mars_tprate_em_i2000-07-01_lt0.tif  
              202407-ECMWF_SEAS-V_i_a-ensemble_mean.tif  
names       : 2000-07-01.0, 2024-07-01.0 
min values  :           ? ,        0.000 
max values  :           ? ,     1741.675

Let’s crop to a smaller extent (Nicaragua) to get a better look at whats going on:

Hard to tell from just plotting raw values

Code 
gdf_nica <-   ne_countries(country = "nicaragua",returnclass = "sf")

r_mars_cropped <- crop(r_mars[[1]],gdf_nica)

plot(r_mars_cropped)

Code 
r_aws_cropped <- crop(r_aws[[1]],gdf_nica)
plot(r_aws_cropped)

Therefore let’s plot the pixel centroids

Raw AWS derived vs MARS

By plotting the pixel centroids we can see that the grid cells of the MARS derived COG are in fact aligned with the original AWS derived COG & properly set CRS AWS derived COG.

Code 
plot(terra::as.points(r_mars_cropped), col = "red")
plot(terra::as.points(r_aws_cropped), col = "blue",add=TRUE, pch=3)

CRS reset AWS derived vs MARS

Code 
r_aws_crs_set_cropped <- crop(r_aws_copy[[1]],gdf_nica)
plot(terra::as.points(r_mars_cropped), col = "red")
plot(terra::as.points(r_aws_crs_set_cropped), col = "blue",add=TRUE, pch=3)

Just as a final gut check add a random raw ensemble from GRIB file to make sure it also lines up. Ensemble 1 is added as green.

Code 
bucket_df <- get_formatted_bucket_df(bucket = Sys.getenv("BUCKET_NAME"))
bucket_df_filt <- bucket_df %>% 
    filter(
        size_mb>12,
        str_detect(date,"2024-07")
        )

key <- bucket_df_filt$Key[1]
filename = bucket_df_filt$filename[1]
td <- tempdir()
tf <- file.path(td, filename)


aws.s3::save_object(bucket=Sys.getenv("BUCKET_NAME"),object = key, file = tf,overwrite = T)
[1] "/var/folders/61/cp06zhcj4y76q7rfx0qlm06c0000gn/T//RtmpD6zqoT/2024-07-05_12-00-47.grib"
Code 
r_grib <- rast(tf)[[1]]

# since we are just adding to the plot we dont actually have to crop
plot(terra::as.points(r_mars_cropped), col = "red")
plot(terra::as.points(r_aws_cropped), col = "blue",add=TRUE, pch=3)
plot(terra::as.points(r_grib), col = "green",add=TRUE, pch=6)

However, the pixel centroids clearly reveal that the AWS derived COG should NOT be reprojected.

Code 
r_aws_proj_cropped <- crop(r_aws_proj[[1]],gdf_nica)
plot(terra::as.points(r_mars_cropped), col = "red")
plot(terra::as.points(r_aws_proj_cropped), col = "blue",add=TRUE, pch=3)

Conclusion

  • The MARS derived COG and the AWS derived COG are in fact aligned and have the same CRS, but the CRS of the AWS derived COG is not properly set as it is somehow getting wrongly imported from the GRIB file.This has not previously been an issue because the AWS CRS was properly set as a post-processing step in the adhoc work,but it would be cleaner to integrate this into pipeline
  • When the AWS NRT Global GRIB -> COG (ensemble mean) pipeline is established it can follow the same CRS setting procedure as the MARS derived COG.