R/aegis_lookup.R

Defines functions aegis_lookup

# TODO:: apply and tapply can be converted into data.table computations for speed
## WARNING:: if your polygons are misspecified, they can result in some au's being dropped .. check your polygons

aegis_lookup = function( 
  pL = NULL, 
  variable_name=NULL,
  LUT = NULL,    # look up table from which to obtain results
  LOCS = NULL,   # look up locations for which results are desired
  LUT_AU=NULL,   # areal units associated with LUT
  LOCS_AU=NULL, 
  project_class="core",   # "project_class" to look up from 
  output_format="points", 
  statvars = c("mean"),
  tz="America/Halifax", 
  FUNC=mean, 
  returntype = "vector",
  space_resolution=NULL,   # spatial planar discretization
  time_resolution=NULL,  # year, seasonal discretization
  year.assessment=NULL, ...
) {


 if (0) {
   # testing using snow data 
    p = bio.snowcrab::snowcrab_parameters(
      project_class="carstm",
      yrs=1999:2021,
      areal_units_type="tesselation",
      carstm_model_label = "default",   # default is the name of areal_units_type
      selection = list(type = "number")
    )

    M = snowcrab.db( p=p, DS="biological_data" )  
   
    M = M[ which( M$yr %in% 2011:2015 ) ,]  # reduce data size
    dim(M)

    sppoly = areal_units( p=p )  # poly to estimate a surface 

    # loadfunctions("aegis")
    # debug(aegis_lookup)
    if (0) {
      # generics
     pL = p
      variable_name=list( "predictions" )
      # LUT = NULL    # look up table from which to obtain results
      LOCS = M[, c("lon", "lat", "timestamp")]   # look up locations for which results are desired
      LUT_AU=NULL   # areal units associated with LUT
      LOCS_AU=sppoly 
      project_class="core"   # "project_class" to look up from 
      output_format="points" 
      statvars = c("mean")
      tz="America/Halifax" 
      FUNC=mean 
      space_resolution=1
    }
    
    loadfunctions("aegis")
    loadfunctions("aegis.survey")
    

    # test raw data lookup
    # spatial
    pL = bathymetry_parameters( project_class="core"  )
    LUT = aegis_survey_lookuptable(
        aegis_project="bathymetry",
        project_class="core",  
        DS="aggregated_data",
        pL=pL 
      )

    # lookup raw data
    o0 = aegis_lookup(  pL=pL, LOCS=M[, c("lon", "lat")],  LUT=LUT,
      project_class="core", output_format="points" , variable_name=c("z.mean" ),
      returntype="sf" 
    ) 
    plot(M[["z"]] ~ o0[["z.mean"]])

    # same as above but return as vector
    o1 = aegis_lookup( pL = pL, LOCS=M[, c("lon", "lat")],  LUT=LUT,
      project_class="core", output_format="points" ,  variable_name=c( "z.mean" ),
      returntype="vector"
    ) 
    plot(M[["z"]] ~ o1)

    # lookup from stmv solution
    pL = bathymetry_parameters( project_class="stmv"  )
    LUT = aegis_survey_lookuptable( aegis_project="bathymetry", project_class="stmv", 
        DS="complete",
        pL=pL )

    o2 = aegis_lookup( pL = pL, LOCS=M[, c("lon", "lat")],  LUT=LUT,  
      project_class="stmv", output_format="points" ,   variable_name=c( "z", "dZ", "ddZ") 
    ) 
    plot(M[["z"]] ~ o2[["z"]])

    # spatial averaging to areal unit
    # this one is slow .. could be sped up if used 
    o3 = aegis_lookup( pL = pL, LOCS=M[, c("lon", "lat")],  LOCS_AU=sppoly, LUT=LUT,  
      project_class="stmv", output_format="areal_units" ,  variable_name=c( "z", "dZ", "ddZ") 
    ) 

    # spatial averaging by areal unit
    o4 = aegis_lookup( pL = pL, LOCS=M[, c("lon", "lat")],  LOCS_AU=sppoly,  LUT=LUT,  
      project_class="core", output_format="areal_units" ,   variable_name=c( "z.mean", "z.sd", "z.n") 
    ) 

    # lookup from carstm areal unit predictions
    pL = bathymetry_parameters( project_class="carstm",  carstm_model_label="default" )
    LUT = aegis_survey_lookuptable(
        aegis_project="bathymetry",
        project_class="carstm",  
        DS="carstm_predictions",
        pL=pL 
      )

    o5 = aegis_lookup( pL = pL, LOCS=M[, c("lon", "lat")],  LUT=LUT,
      project_class="carstm", output_format="points", variable_name=list( "predictions" ), statvars=c("mean", "sd"), space_resolution=min(p$pres)
    ) 
    plot(M[["z"]] ~ o5)

    if (0){
    pL = bathymetry_parameters( project_class="carstm",  carstm_model_label="default" )
      # to get spatial random effects
        LUT = aegis_survey_lookuptable(
            aegis_project="bathymetry",
            project_class="carstm",  
            DS="carstm_randomeffects",
            pL=pL 
          )

        o5r = aegis_lookup( pL = pL, LOCS=M[, c("lon", "lat")],  LUT=LUT,
          project_class="carstm", output_format="points", variable_name=list( c( "space", "re_total") ), statvars=c("mean", "sd"), space_resolution=min(p$pres)
        ) 
    }
    # slow 
    o6 = aegis_lookup( pL = pL, LOCS=M[, c("lon", "lat")],  LOCS_AU=sppoly, LUT=LUT,
      project_class="carstm", output_format="areal_units", variable_name=list( "predictions"  ), statvars=c("mean", "sd"), space_resolution=min(p$pres) /2
    ) 

    # au to au match (LOCS=NULL)
    o7 = aegis_lookup( pL = pL, LOCS=list(AUID=sppoly$AUID), LOCS_AU=sppoly,  LUT=LUT,
      project_class="carstm", output_format="areal_units", variable_name=list( "predictions"  ), statvars=c("mean", "sd"), space_resolution=min(p$pres) /2
    ) 

    # au to au match (LOCS=NULL)
    o7 = aegis_lookup( pL = pL, LOCS=sppoly[,"AUID"], LOCS_AU=sppoly, LUT=LUT,
      project_class="carstm", output_format="areal_units", variable_name=list( "predictions"  ), statvars=c("mean", "sd"), space_resolution=min(p$pres) /2
    ) 

    # au to au match (LOCS=NULL)
    o7 = aegis_lookup( pL = pL, LOCS=sppoly$AUID, LOCS_AU=sppoly, LUT=LUT,
      project_class="carstm", output_format="areal_units", variable_name=list( "predictions"  ), statvars=c("mean", "sd"), space_resolution=min(p$pres) /2
    ) 

    # consistency checks
    plot(o1 ~o2$z)
    plot(o1 ~o3$z)
    plot(o1 ~o4$z)
    plot(o1 ~o5$predictions_mean) 
    plot(o1 ~o6$predictions_mean) 


    # space-time
    year.assessment=2023
    pL = speciescomposition_parameters( 
        project_class="carstm", yrs=1999:year.assessment, carstm_model_label="default" ,
        variabletomodel="pca1"
    )
     
        LUT = aegis_survey_lookuptable(
            aegis_project="speciescomposition",
            project_class="carstm",  
            DS="carstm_predictions",
            pL=pL 
          )

 
    # o1 = aegis_lookup(  pL=pL, LOCS=M[, c("lon", "lat", "timestamp")],   LUT=LUT,
    #   project_class="core", output_format="points" , variable_name=c( "pca1"), space_resolution=2  
    # ) 

    # o1 = aegis_lookup(  pL=pL, LOCS=M[, c("lon", "lat", "timestamp")],     LUT=LUT,
    #   project_class="core", output_format="points" , variable_name=c( "pca1"), space_resolution=2  
    # ) 

    # o2 = aegis_lookup(  pL=pL, LOCS=M[, c("lon", "lat", "timestamp")],     LUT=LUT,
    #   project_class="stmv", output_format="points" ,   variable_name=c( "pca1") 
    # ) 

    # o3 = aegis_lookup(  pL=pL, LOCS=M[, c("lon", "lat", "timestamp")], LOCS_AU=sppoly,   LUT=LUT,
    #   project_class="core", output_format="areal_units" ,  variable_name=list( "pca1",  "pca2", "ca1", "ca2" )
    # ) 
 
    o4 = aegis_lookup(  pL=pL, LOCS=M[, c("lon", "lat", "timestamp")],   LUT=LUT,
      project_class="carstm", output_format="points", variable_name=list( "predictions" ), statvars=c("mean", "sd")
    ) 

    o5 = aegis_lookup(  pL=pL, LOCS=M[, c("lon", "lat", "timestamp")], LOCS_AU=sppoly,    LUT=LUT,
      project_class="carstm", output_format="areal_units" , variable_name=list( "predictions" ), statvars=c("mean", "sd"), space_resolution=min(p$pres)/2
    ) 

    mm = expand.grid(AUID=sppoly$AUID, timestamp=lubridate::date_decimal( p$yrs, tz="America/Halifax" ))
    mm = expand.grid(AUID=sppoly$AUID, timestamp= p$yrs )

    o6 = aegis_lookup(  pL=pL, LOCS=mm, LOCS_AU=sppoly,    LUT=LUT,
      project_class="carstm", output_format="areal_units" , variable_name=list( "predictions" ), statvars=c("mean", "sd"), space_resolution=min(p$pres)/2
    ) 


    # consistency checks
    plot(o1$pca1.mean ~o3$pca1.mean)
    plot(o1$pca1.mean ~o4$predictions_mean) 
    plot(o1$pca1.mean ~o5$predictions_mean) 


    # space-time-season

    year.assessment=2023
    pL = temperature_parameters( 
        project_class="carstm", yrs=1999:year.assessment, carstm_model_label="default" 
    )
     LUT = aegis_survey_lookuptable(
            aegis_project="temperature",
            project_class="carstm",  
            DS="carstm_predictions",
            pL=pL 
          )

 
    # o1 = aegis_lookup(  pL=pL, LOCS=M[, c("lon", "lat", "timestamp")],    LUT=LUT, 
    #   project_class="core", output_format="points" ,  variable_name=c( "t.mean", "t.sd", "t.n"), space_resolution=2  
    # ) 

    # o2 = aegis_lookup(  pL=pL, LOCS=M[, c("lon", "lat", "timestamp")],    LUT=LUT,  
    #   project_class="stmv", output_format="points" ,  variable_name=c( "t") 
    # ) 

    # o3 = aegis_lookup(  pL=pL, LOCS=M[, c("lon", "lat", "timestamp")], LOCS_AU=sppoly,    LUT=LUT,
    #   project_class="core", output_format="areal_units" ,  variable_name=list( "t.mean",  "t.sd", "t.n"  )
    # ) 
 
    o4 = aegis_lookup(  pL=pL, LOCS=M[, c("lon", "lat", "timestamp")], LUT=LUT,
      project_class="carstm", output_format="points", variable_name=list( "predictions" ), statvars=c("mean", "sd")
    ) 
    plot(M$t ~o4 )

    o5 = aegis_lookup(  pL=pL, LOCS=M[, c("lon", "lat", "timestamp")], LOCS_AU=sppoly,  LUT=LUT,
      project_class="carstm", output_format="areal_units" , variable_name=list( "predictions" ), statvars=c("mean", "sd"), space_resolution=min(p$pres)/2
    ) 


   # mm = expand.grid(AUID=sppoly$AUID, timestamp=lubridate::date_decimal( p$yrs, tz="America/Halifax" ))
    mm = expand.grid(AUID=sppoly$AUID, timestamp= p$yrs )

    o6 = aegis_lookup(  pL=pL, LOCS=mm, LOCS_AU=sppoly,  LUT=LUT,
      project_class="carstm", output_format="areal_units" , variable_name=list( "predictions" ), statvars=c("mean", "sd"), space_resolution=min(p$pres)/2
    ) 

    # consistency checks
    plot(o1$t.mean ~o3$t.mean)
    plot(o1$t.mean ~o4$predictions_mean) 
    plot(o1$t.mean ~o5$predictions_mean) 



    if (0) {
      p = bathymetry_parameters(  project_class="stmv"  )
      LUT = bathymetry_db ( p=p, DS="complete" )   
    }
 
  }  ## END consistency checks

  require(stars)

  # -----------
  if (is.null(LUT)) stop("LUT is required")
  
  i = attributes(LUT)
  
  if (is.null(pL)) {
     if ("pL" %in% names(i)) pL = i[["pL"]]
  }

  if (is.null(pL)) stop("pL is required")


  if (is.null(space_resolution)) {
      if (exists( "pres", pL)) {
        space_resolution = pL$pres
      }
  }
  
  if (is.null(time_resolution)) {
      if (exists( "tres", pL)) {
        time_resolution = pL$tres
      }
  }
      
  if (is.null(space_resolution)) {
    if (length(i)>0) {
      if ("space_resolution" %in% names(i)) space_resolution = i[["space_resolution"]]
    }
  }

  if (is.null(time_resolution)) { 
    if (length(i)>0) {
      if ("time_resolution" %in% names(i)) time_resolution = i[["time_resolution"]]
    }
  }

  if (is.vector(LOCS)) LOCS = list(AUID=LOCS)
  if (!is.null(LOCS)) setDT(LOCS)

  crs_lonlat =  st_crs(projection_proj4string("lonlat_wgs84"))
  

  if (is.null(variable_name)) {
    if ( project_class %in% c("carstm" )) {
      variable_name = LUT$fm$dependent_variable
    } else {
      variable_name = setdiff( names(LUT), c("plon", "plat", "lon", "lat", "timestamp", "year", "dyear", "yr" ))
    }
  } 


  if ( project_class %in% c("core", "stmv", "hybrid") )  {
    
    if ( pL$dimensionality == "space" ) {
        if ( space_resolution != pL$pres ) {
          setDT(LUT)
          # regrid to another resolution
          LUT$plon = trunc(LUT$plon / space_resolution + 1 ) * space_resolution
          LUT$plat = trunc(LUT$plat / space_resolution + 1 ) * space_resolution
          LUT = LUT[, setNames(.(mean( get(variable_name), na.rm=TRUE) ), variable_name), by=list(plon, plat) ]
          setDF(LUT)
        }
        attr(LUT, "space_resolution") = space_resolution
        attr(LUT, "pL") = pL
    }

  
    if (pL$dimensionality == "space-time"  ) {
    
        if ( time_resolution != pL$tres ) {
          setDT(LUT)
          LUT$yr = trunc( (LUT$yr + LUT$dyear )/ time_resolution + 1 ) * time_resolution
        }
        if ( space_resolution != pL$pres ) {
          # regrid to another resolution
          setDT(LUT)
          LUT$plon = trunc(LUT$plon / space_resolution + 1 ) * space_resolution
          LUT$plat = trunc(LUT$plat / space_resolution + 1 ) * space_resolution
        }
        if ( time_resolution != pL$tres | space_resolution != pL$pres ) {
          setDT(LUT)
          LUT = LUT[, setNames(.(mean( get(variable_name), na.rm=TRUE) ), variable_name), by=list(plon, plat, yr ) ]
          LUT$timestamp = lubridate::date_decimal( LUT$yr, tz=tz )
        }
        setDF(LUT)
        attr(LUT, "space_resolution") = space_resolution
        attr(LUT, "time_resolution") = time_resolution
        attr(LUT, "pL") = pL
    }

    if (pL$dimensionality == "space-time-cyclic")  {
    
        if ( time_resolution != pL$tres ) {
          setDT(LUT)
          LUT$dyear = trunc(LUT$dyear / time_resolution + 1 ) * time_resolution
        }
        if ( space_resolution != pL$pres ) {
          # regrid to another resolution
          setDT(LUT)
          LUT$plon = trunc(LUT$plon / space_resolution + 1 ) * space_resolution
          LUT$plat = trunc(LUT$plat / space_resolution + 1 ) * space_resolution
        }
        if ( time_resolution != pL$tres | space_resolution != pL$pres ) {
          setDT(LUT)
          LUT = LUT[, setNames(.(mean( get(variable_name), na.rm=TRUE) ), variable_name), by=list(plon, plat, yr, dyear) ]
          LUT$timestamp = lubridate::date_decimal( LUT$yr+LUT$dyear, tz=tz )
        }
        setDF(LUT)
        attr(LUT, "space_resolution") = space_resolution
        attr(LUT, "time_resolution") = time_resolution
        attr(LUT, "pL") = pL
      }

  }


  # ------------------

  if (pL$dimensionality =="space" ) {

    if ( project_class %in% c("core", "stmv", "hybrid") & output_format == "points" )  {
      
      if (!exists("plon", LUT))  LUT = lonlat2planar(LUT, proj.type=pL$aegis_proj4string_planar_km)
      if (!exists("plon", LOCS)) LOCS = lonlat2planar(LOCS, proj.type=pL$aegis_proj4string_planar_km) # get planar projections of lon/lat in km

      LOCS$tplon = trunc(LOCS$plon / space_resolution + 1 ) * space_resolution
      LOCS$tplat = trunc(LOCS$plat / space_resolution + 1 ) * space_resolution

      plon_range = range( LUT$plon )
      plat_range = range( LUT$plat )
      plons = seq(plon_range[1], plon_range[2], by=space_resolution)
      plats = seq(plat_range[1], plat_range[2], by=space_resolution)
      Sdims = c(length(plons), length(plats))
      Sres = c(space_resolution, space_resolution)
      Sorigin = c( plon_range[1], plat_range[1] )

      ii = match( 
          array_map( "xy->1", LOCS[, c("tplon","tplat")], dims=Sdims, res=Sres, origin=Sorigin ),
          array_map( "xy->1", LUT[,c("plon","plat")], dims=Sdims, res=Sres, origin=Sorigin )
      )

      LOCS$tplon = LOCS$tplat  = NULL
      gc()

      for (vnm in variable_name) {
        if ( vnm %in% names(LUT )) {
          LOCS[[vnm]] = NA  
          LOCS[[vnm]] = LUT[ ii, vnm ]
        } 
      } 

    }

    if ( project_class %in% c("core", "stmv", "hybrid") & output_format == "areal_units" )  {

      LOCS0 = NULL
      if (exists("AUID", LOCS)) {
        # this is a special class of lookup taking LUT (usually empirical or aggregated data) and estimating from it for an areal unit (via rasterization)
        LOCS0 = LOCS # store as modifications will be made to LOCS

        # next, get a x-y coordinate for each LOC by rasterizing to LOCS_AU  
        LOCS_AU = st_make_valid(LOCS_AU)
        LOCS_AU = sf::st_transform( LOCS_AU, crs=st_crs(pL$aegis_proj4string_planar_km) )

        LOCS_AU$au_index = 1:nrow(LOCS_AU)
        LOCS_AU_raster = stars::st_rasterize( LOCS_AU["au_index"], dx=space_resolution, dy=space_resolution )

        # overwrite LOCS with discretized points 
        LOCS = sf::st_as_sf( LOCS_AU_raster, as_points=TRUE, na.rm=FALSE )
        st_crs(LOCS) = st_crs( LOCS_AU )
        LOCS$AUID = LOCS_AU$AUID[ st_extract(LOCS_AU_raster, LOCS )$au_index ]

      } else {

        if (exists("lon", LOCS)) {
          LOCS = sf::st_as_sf( LOCS, coords=c("lon", "lat") )
          st_crs(LOCS) = st_crs( projection_proj4string("lonlat_wgs84") )
        } else if (exists("plon", LOCS)) {
          LOCS = sf::st_as_sf( LOCS, coords=c("plon", "plat") )
          st_crs(LOCS) = st_crs(pL$aegis_proj4string_planar_km) 
        } 

      }

      if (exists("lon", LUT)) {
        LUT = sf::st_as_sf( LUT, coords=c("lon", "lat") )
        st_crs(LUT) = st_crs( projection_proj4string("lonlat_wgs84") )
      } else if (exists("plon", LUT)) {
        LUT = sf::st_as_sf( LUT, coords=c("plon", "plat") )
        st_crs(LUT) = st_crs(pL$aegis_proj4string_planar_km) 
      } 

      if (!exists("AUID", LOCS_AU))  {
        message ("AUID not found in the polygons, setting AUID as row number of polygons")
        LOCS_AU$AUID = as.character(1:nrow(LOCS_AU))
      }
      LOCS_AU = sf::st_transform( LOCS_AU, crs=st_crs(LUT) )  # 


      variable_name = intersect( names(LUT), variable_name )
      LUT = LUT[, variable_name]
      LU_map =  st_points_in_polygons( pts=LUT, polys=LOCS_AU[, "AUID"], varname= "AUID" ) 

      if (!exists("AUID", LOCS ))  LOCS[["AUID"]]  = st_points_in_polygons( pts=LOCS, polys=LOCS_AU[, "AUID"], varname= "AUID" ) 

      st_geometry(LUT) = NULL
      setDT(LUT)
      for (vnm in variable_name) {
        if ( vnm %in% names(LUT )) {  
          # LUT_regridded = tapply( st_drop_geometry(LUT)[, vnm], LU_map, FUN=FUNC, na.rm=TRUE )
          LUT_regridded = LUT[, setNames(.(mean( get(vnm), na.rm=TRUE) ), vnm), by=.(LU_map) ]
          LOCS[[ vnm ]] = LUT_regridded[ match( LOCS$AUID, LUT_regridded$LU_map ), vnm, with=FALSE ]
        }   
      }

      # revert LOCS to original structure and add the vars to it
      if (!is.null(LOCS0)) {
        # summarize LOCS to to original AUID's
        if ( "sf" %in% class(LOCS) ) st_geometry(LOCS) = NULL
        LOCS = as.data.table( LOCS )
        # o = LOCS[,  lapply(.SD, mean, na.rm=TRUE), by=AUID, .SDcols=variable_name]   
        o = LOCS[,  setNames(.(mean( get(variable_name), na.rm=TRUE) ), variable_name), by=AUID ]   
        LOCS0[,variable_name] = o[ match( LOCS0$AUID, o$AUID), .SD, .SDcols=variable_name  ]        
        LOCS = LOCS0
      }

    }


    if ( project_class %in% c("carstm" ) & output_format == "points" )  {

      if (is.null(LUT_AU)) LUT_AU = LUT$sppoly
      LUT_AU = st_make_valid(LUT_AU)
      LUT_AU = st_cast( LUT_AU, "MULTIPOLYGON")
      # LUT_AU = st_cast(LUT_AU, "POLYGON" )
      LUT_AU = st_make_valid(LUT_AU)
      LUT_AU = sf::st_transform( LUT_AU, crs=st_crs(pL$aegis_proj4string_planar_km) )

      if (!exists("AUID", LUT_AU)) {
        message ("AUID not found in LOCS_AU polygons, setting AUID as row number of polygons")
        LUT_AU$AUID = as.character(LUT_AU$au_index)
      }
      bm = match( LUT_AU$AUID, LUT$space_name )  
            
      if (exists("lon", LOCS)) {
        LOCS = sf::st_as_sf( LOCS, coords=c("lon", "lat") )
        st_crs(LOCS) = st_crs( projection_proj4string("lonlat_wgs84") )
        LOCS = sf::st_transform( LOCS, crs=st_crs(pL$aegis_proj4string_planar_km) )
      } else if  (exists("plon", LOCS)) { 
        LOCS = sf::st_as_sf( LOCS, coords=c("lon", "lat") )
        st_crs(LOCS) = st_crs(pL$aegis_proj4string_planar_km)
      }

      LOCS$AUID = st_points_in_polygons( pts=LOCS, polys = LUT_AU[, "AUID"], varname= "AUID" )   
      ii = match( LOCS$AUID, LUT$space_name[bm] ) 

      for (g in 1:length(statvars)) {
        stat_var = statvars[g]
        for (h in 1:length(variable_name) ) {
          vnh = variable_name[[h]] 
          if (exists(vnh[[1]], LUT)) {
            vnm = paste( paste0(vnh,  collapse="_"),  stat_var, sep="_" )
            LUT_AU[[vnm]] = carstm_results_unpack( LUT, vnh ) [ bm, stat_var ]
            LOCS[[vnm]] = LUT_AU[[vnm]] [ ii ]
          }
        }
      } 

    }



    if ( project_class %in% c("carstm") & output_format == "areal_units" )  {
      
      if (is.null(LUT_AU)) LUT_AU = LUT$sppoly
      LUT_AU = st_make_valid(LUT_AU)
      LUT_AU = st_cast( LUT_AU, "MULTIPOLYGON")
      # LUT_AU = st_cast(LUT_AU, "POLYGON" )
      LUT_AU = st_make_valid(LUT_AU)
      LUT_AU = sf::st_transform( LUT_AU, crs=st_crs(pL$aegis_proj4string_planar_km) )

      if (!exists("AUID", LUT_AU)) LUT_AU$AUID = as.character(1:nrow(LUT_AU))
      bm = match( LUT_AU$AUID, LUT$space_name )  # should not be required but just in case things are subsetted

      # lut_uid is A LOCAL index of LUT_AU /LUT .. rasterize it
      LUT_AU$lut_uid = 1:nrow(LUT_AU)
      LUT_AU_raster = stars::st_rasterize( LUT_AU["lut_uid"], dx=space_resolution, dy=space_resolution )
      LUT_AU_pts = sf::st_as_sf( LUT_AU_raster, as_points=TRUE, na.rm=FALSE )
      st_crs(LUT_AU_pts) = st_crs( LUT_AU ) 

      # format output polygon structure 
      LOCS_AU = sf::st_transform( LOCS_AU, crs=st_crs(LUT_AU) )  # LOCS_AU .... must be sent ... <---------
      LOCS_AU = st_make_valid(LOCS_AU)
      LOCS_AU = st_buffer( LOCS_AU, 0 )
      LOCS_AU = st_cast( LOCS_AU, "MULTIPOLYGON")
#        LOCS_AU = st_cast( LOCS_AU, "POLYGON")
      if (!exists("AUID", LOCS_AU))  {
        message ("AUID not found in the polygons, setting AUID as row number of polygons")
        LOCS_AU$AUID = as.character(1:nrow(LOCS_AU))
      }

      if (! inherits(LOCS, "sf") )   {
        if (exists("lon", LOCS )) {
          LOCS = st_as_sf( LOCS, coords=c("lon","lat"), crs=st_crs(projection_proj4string("lonlat_wgs84")) )
          LOCS = sf::st_transform( LOCS, crs=st_crs(pL$aegis_proj4string_planar_km) )
          if (!exists("AUID", LOCS )) LOCS[["AUID"]]  = st_points_in_polygons( pts=LOCS, polys=LOCS_AU[, "AUID"], varname= "AUID" ) 
        }
      }

      # covert LUT_AU to LOCS_AU
      # id membership of LOCS in LOCS_AU

      # id membership of LUT raster in LOCS_AU
      LUT_AU_pts_LOCS_AU_AUID = st_points_in_polygons( pts=LUT_AU_pts, polys=LOCS_AU[,"AUID"], varname="AUID" ) 

      pts_AU = match(LUT_AU_pts$lut_uid, LUT_AU$lut_uid[match(LUT$space_name, LUT_AU$AUID)] ) ## (layer==lut_uid) -- to -- LUT

      aggFUN = function( LUV ) tapply(X=LUV, INDEX=LUT_AU_pts_LOCS_AU_AUID, FUN=FUNC, na.rm=TRUE) 

      for (g in 1:length(statvars)) {
        stat_var = statvars[g]
        for (h in 1:length(variable_name) ) {
          vnh = variable_name[[h]] 
          if (exists(vnh[[1]], LUT)) {
            vnm = paste( paste0(vnh,  collapse="_"),  stat_var, sep="_" )
            o = carstm_results_unpack( LUT, vnh ) 
            kk = length(dim(o))
            if (kk == 2) {
              ov = o[ pts_AU,  which(dimnames(o)$stat == stat_var)  ] 
              LUT_as_LOCS_AU = aggFUN( ov)
              space_index = match( as.character(LOCS$AUID), as.character(dimnames(LUT_as_LOCS_AU )[[1]] )  )   # AUID of LOCS_AU (from LUT_AU_pts_AUID)
              if (all(is.na(space_index))) {
                LOCS[,vnm] = LUT_as_LOCS_AU  # fiddling required when only 1 AU
              } else {
                LOCS[,vnm] = LUT_as_LOCS_AU[ space_index  ]
              }
            }
          }
        } 
      }
    }

  }


  ######################################################

  if (pL$dimensionality =="space-time" ) {


    if ( project_class %in% c("core", "stmv", "hybrid") & output_format == "points" )  {

      if (!exists("plon", LUT))  LUT = lonlat2planar(LUT, proj.type=pL$aegis_proj4string_planar_km)
      if (!exists("plon", LOCS)) LOCS = lonlat2planar(LOCS, proj.type=pL$aegis_proj4string_planar_km) # get planar projections of lon/lat in km
      if (!exists("dyear", LOCS) | (!exists("yr", LOCS)) ) {
        LOCS$tiyr = lubridate::decimal_date( LOCS$timestamp  ) 
        LOCS$yr = trunc( LOCS$tiyr )
        LOCS$dyear = LOCS$tiyr - LOCS$yr   
      } 

      LOCS$tplon = trunc(LOCS$plon / space_resolution + 1 ) * space_resolution
      LOCS$tplat = trunc(LOCS$plat / space_resolution + 1 ) * space_resolution
      LOCS$tdyear = trunc(LOCS$dyear / time_resolution + 1 ) * time_resolution

      nY = diff( range(LUT$yr) )
      nW = round(1/time_resolution)
      Tdims = c( nY, nW )
      Tres = c(1, time_resolution)
      Torigin = c( min(LUT$yr), 0) 

      plon_range = range( LUT$plon )
      plat_range = range( LUT$plat )
      plons = seq(plon_range[1], plon_range[2], by=space_resolution)
      plats = seq(plat_range[1], plat_range[2], by=space_resolution)
      Sdims = c(length(plons), length(plats))
      Sres = c(space_resolution, space_resolution)
      Sorigin = c( plon_range[1], plat_range[1] )

      ii = match( 
        paste(
          array_map( "xy->1", LOCS[, c("tplon","tplat")], dims=Sdims, res=Sres, origin=Sorigin ), 
          array_map( "ts->1", LOCS[, c("yr", "tdyear")], dims=Tdims, res=Tres, origin=Torigin ), 
          sep="_"
        ), 
        paste( 
          array_map( "xy->1", LUT[,c("plon","plat")], dims=Sdims, res=Sres, origin=Sorigin ), 
          array_map( "ts->1", LUT[,c("yr", "dyear")], dims=Tdims, res=Tres, origin=Torigin ), 
          sep="_" 
        ) 
      )

      LOCS$tplon = LOCS$tplat = LOCS$tdyear = NULL
      gc()

      for (vnm in variable_name) {
        if ( vnm %in% names(LUT )) {
          LOCS[[vnm]] = NA
          LOCS[[vnm]] = LUT[ ii, vnm ]

        } 

      }
  
    }


    if ( project_class %in% c("core", "stmv", "hybrid") & output_format == "areal_units" )  {

      if (exists("lon", LUT)) {
        LUT = sf::st_as_sf( LUT, coords=c("lon", "lat") )
        st_crs(LUT) = st_crs( projection_proj4string("lonlat_wgs84") )
      } else if (exists("plon", LUT)) {
        LUT = sf::st_as_sf( LUT, coords=c("plon", "plat") )
        st_crs(LUT) = st_crs(pL$aegis_proj4string_planar_km) 
      } 

      Tdims = c( diff( range(LUT$yr) ), pL$nw )
      Tres = c(1, pL$tres)
      Torigin = c( min(LUT$yr), 0 ) 

      Sdims = c( pL$nplons, pL$nplats )
      Sres = c( pL$pres, pL$pres )
      Sorigin = pL$origin

      if (!exists("AUID", LOCS_AU))  {
        message ("AUID not found in LOCS_AU polygons, setting AUID as row number of polygons")
        LOCS_AU$AUID = as.character(1:nrow(LOCS_AU))
      }

      if ( ! "sf" %in% class(LOCS) ) {
        LOCS = st_as_sf( LOCS, coords=c("lon","lat"), crs=st_crs(projection_proj4string("lonlat_wgs84")) )
      }

      LOCS_AU = sf::st_transform( LOCS_AU, crs=st_crs(LUT) )

      if (! inherits(LOCS$timestamp, "POSIXct") )  LOCS$timestamp =  lubridate::date_decimal( LOCS$timestamp, tz=tz )
      if (! exists("yr", LOCS) ) LOCS$yr = lubridate::year(LOCS$timestamp) 

      LU_map = paste( 
        st_points_in_polygons( pts=LUT, polys=LOCS_AU[, "AUID"], varname= "AUID" ), 
        array_map( "ts->1", LUT[,c("yr", "dyear")], dims=Tdims, res=Tres, origin=Torigin ), 
        sep="_"
      )

      LOCS_map =  paste( 
        st_points_in_polygons( pts=LOCS, polys = LOCS_AU[, "AUID"], varname= "AUID" ),  
        array_map( "ts->1", LOCS[, c("yr", "tdyear")], dims=Tdims, res=Tres, origin=Torigin ), 
        sep="_"
      )
  
      LUT = setDT(st_drop_geometry(LUT))

      for (vnm in variable_name) {
        if ( vnm %in% names(LUT )) {  
          LUT_regridded = LUT[, setNames(.(mean( get(vnm), na.rm=TRUE) ), vnm), by=.(LU_map) ]
          LOCS[[ vnm ]] = LUT_regridded[ match( LOCS$AUID, LUT_regridded$LU_map ), vnm, with=FALSE ]
        }   
      }


    }




    if ( project_class %in% c("carstm" ) & output_format == "points" )  {
      
      # discretize LUT_AU to LOCS_AU and then re-estimate means by LOCS_AU$AUID then tag each estimate to AUID of LOCS

      ny = length(LUT$time_name)
      yr0 = min(as.numeric(LUT$time_name))

      if (is.null(LUT_AU)) LUT_AU = LUT$sppoly
      LUT_AU = st_make_valid(LUT_AU)
      LUT_AU = st_cast( LUT_AU, "MULTIPOLYGON")
      # LUT_AU = st_cast(LUT_AU, "POLYGON" )
      LUT_AU = st_make_valid(LUT_AU)
      LUT_AU = sf::st_transform( LUT_AU, crs=st_crs(pL$aegis_proj4string_planar_km) )
      LUT_AU$au_index = 1:nrow(LUT_AU)
      if (!exists("AUID", LUT_AU)) {
        message ("AUID not found in LOCS_AU polygons, setting AUID as row number of polygons")
        LUT_AU$AUID = as.character(LUT_AU$au_index)
      }
      bm = match( LUT_AU$AUID, LUT$space_name )  
            
      if (exists("lon", LOCS)) {
        LOCS = sf::st_as_sf( LOCS, coords=c("lon", "lat") )
        st_crs(LOCS) = st_crs( projection_proj4string("lonlat_wgs84") )
        LOCS = sf::st_transform( LOCS, crs=st_crs(pL$aegis_proj4string_planar_km) )
      } else if  (exists("plon", LOCS)) { 
        LOCS = sf::st_as_sf( LOCS, coords=c("lon", "lat") )
        st_crs(LOCS) = st_crs(pL$aegis_proj4string_planar_km)
      }

      LOCS$AUID = st_points_in_polygons( pts=LOCS, polys = LUT_AU[, "AUID"], varname= "AUID" )   
      
      if (! inherits(LOCS$timestamp, "POSIXct") )  LOCS$timestamp =  lubridate::date_decimal( LOCS$timestamp, tz=tz )
      if (! exists("yr", LOCS) ) LOCS$yr = lubridate::year(LOCS$timestamp) 
    
    
      ii = cbind( match( LOCS$AUID, LUT$space_name[bm] ) )

      jj = cbind( 
        ii, 
        array_map( "ts->year_index", LOCS[["yr"]], dims=c(ny ), res=c( 1  ), origin=c( yr0 ) ) 
      )

      for (g in 1:length(statvars)) {
        stat_var = statvars[g]
        for (h in 1:length(variable_name) ) {
          vnh = variable_name[[h]] 
          if (exists(vnh[[1]], LUT)) {
            vnm = paste( paste0(vnh,  collapse="_"),  stat_var, sep="_" )
            o = carstm_results_unpack( LUT, vnh ) 
            kk = length(dim(o))
            if (kk == 2) LOCS[[vnm]] = o[ cbind(ii, which(dimnames(o)$stat == stat_var) )  ] 
            if (kk == 3) LOCS[[vnm]] = o[ cbind(jj, which(dimnames(o)$stat == stat_var) )  ] 
          }
        } 
      }

    }


    if ( project_class %in% c("carstm") & output_format == "areal_units" )  {

      # discretize LUT_AU to LOCS_AU and then re-estimate means by LOCS_AU$AUID  then tag each estimate to AUID of LOCS

      # from source data: LUT = modelled predictions; LUT_AU are associated areal units linked by "AUID" 
      #    nw = length(LUT$dyear)

      # format coordinate systems   
      if (is.null(LUT_AU)) LUT_AU = LUT$sppoly
      LUT_AU = st_make_valid(LUT_AU)
      LUT_AU = st_cast( LUT_AU, "MULTIPOLYGON")
      # LUT_AU = st_cast(LUT_AU, "POLYGON" )
      LUT_AU = st_make_valid(LUT_AU)
      LUT_AU = sf::st_transform( LUT_AU, crs=st_crs(pL$aegis_proj4string_planar_km) )
      
      # dims
      ny = length(LUT$time_name)
      yr0 = min(as.numeric(LUT$time_name))

      if (!exists("AUID", LUT_AU)) LUT_AU$AUID = as.character(1:nrow(LUT_AU))
      bm = match( LUT_AU$AUID, LUT$space_name )  # should not be required but just in case things are subsetted

      # lut_uid is A LOCAL index of LUT_AU /LUT .. rasterize it
      LUT_AU$lut_uid = 1:nrow(LUT_AU)
      LUT_AU_raster = stars::st_rasterize( LUT_AU["lut_uid"], dx=space_resolution, dy=space_resolution )
      LUT_AU_pts = sf::st_as_sf( LUT_AU_raster, as_points=TRUE, na.rm=FALSE )
      st_crs(LUT_AU_pts) = st_crs( LUT_AU ) 

      # format output polygon structure 
      LOCS_AU = sf::st_transform( LOCS_AU, crs=st_crs(LUT_AU) )  # LOCS_AU .... must be sent ... <---------
      LOCS_AU = st_make_valid(LOCS_AU)
      LOCS_AU = st_buffer( LOCS_AU, 0 )
      LOCS_AU = st_cast( LOCS_AU, "MULTIPOLYGON")
#        LOCS_AU = st_cast( LOCS_AU, "POLYGON")
      if (!exists("AUID", LOCS_AU))  {
        message ("AUID not found in the LOCS_AU polygons, setting AUID as row number of polygons")
        LOCS_AU$AUID = as.character(1:nrow(LOCS_AU))
      }

      if (! inherits(LOCS, "sf") )   {
        if (exists("lon", LOCS )) {
          LOCS = st_as_sf( LOCS, coords=c("lon","lat"), crs=st_crs(projection_proj4string("lonlat_wgs84")) )
          LOCS = sf::st_transform( LOCS, crs=st_crs(pL$aegis_proj4string_planar_km) )
          if (!exists("AUID", LOCS )) LOCS[["AUID"]]  = st_points_in_polygons( pts=LOCS, polys=LOCS_AU[, "AUID"], varname= "AUID" ) 
        }
      }

      if (! inherits(LOCS$timestamp, "POSIXct") )  LOCS$timestamp =  lubridate::date_decimal( LOCS$timestamp, tz=tz )
      if (! exists("yr", LOCS) ) LOCS$yr = lubridate::year(LOCS$timestamp) 
      
      TIMESTAMP_index1 = array_map( "ts->year_index", LOCS[["yr"]], dims=c(ny ), res=c( 1  ), origin=c( yr0 ) )
      if (any( TIMESTAMP_index1 < 0)) {
        TIMESTAMP_index1[ which( TIMESTAMP_index1 <= 0 ) ] = NA
      }

      # id membership of LUT raster in LOCS_AU
      LUT_AU_pts_LOCS_AU_AUID = st_points_in_polygons( pts=LUT_AU_pts, polys=LOCS_AU[,"AUID"], varname="AUID" ) 

      pts_AU = match(LUT_AU_pts$lut_uid, LUT_AU$lut_uid[match(LUT$space_name, LUT_AU$AUID)] ) ## (layer==lut_uid) -- to -- LUT
      aggFUN = function( LUV ) {
          tapply(X=LUV, INDEX=LUT_AU_pts_LOCS_AU_AUID, FUN=FUNC, na.rm=TRUE) 
      }

      for (g in 1:length(statvars)) {
        stat_var = statvars[g]
        for (h in 1:length(variable_name) ) {
          vnh = variable_name[[h]] 
          if (exists(vnh[[1]], LUT)) {
            vnm = paste( paste0(vnh,  collapse="_"),  stat_var, sep="_" )
            o = carstm_results_unpack( LUT, vnh ) 
            kk = length(dim(o))
            if (kk == 2) {
              ov = o[ pts_AU,  which(dimnames(o)$stat == stat_var)  ] 
              LUT_as_LOCS_AU = aggFUN( ov)
              space_index = match( as.character(LOCS$AUID), as.character(dimnames(LUT_as_LOCS_AU )[[1]] )  )   # AUID of LOCS_AU (from LUT_AU_pts_AUID)
              if (all(is.na(space_index))) {
                LOCS[,vnm] = LUT_as_LOCS_AU  # fiddling required when only 1 AU
              } else {
                LOCS[,vnm] = LUT_as_LOCS_AU[  space_index  ]
              }
            }
            if (kk == 3) {
              ov = o[ pts_AU,, which(dimnames(o)$stat == stat_var)  ] 
              LUT_as_LOCS_AU = apply( ov, MARGIN=c(2), FUN=aggFUN )
              space_index = match( as.character(LOCS$AUID), as.character(dimnames(LUT_as_LOCS_AU )[[1]] )  )   # AUID of LOCS_AU (from LUT_AU_pts_AUID)
              if (all(is.na(space_index))) {
                LOCS[,vnm] = LUT_as_LOCS_AU[ TIMESTAMP_index1  ]  # fiddling required when only 1 AU
              } else {
                LOCS[,vnm] = LUT_as_LOCS_AU[ cbind( space_index, TIMESTAMP_index1 ) ]
              }
            }
          }
        } 
      }
    }

  }  # end dimension
  
  ######################################################


  if (pL$dimensionality =="space-time-cyclic" ) {

    if ( project_class %in% c("core", "stmv", "hybrid") & output_format == "points" )  {

      if (!exists("plon", LOCS))  LOCS = lonlat2planar(LOCS, proj.type=pL$aegis_proj4string_planar_km) # get planar projections of lon/lat in km
      if (!exists("plon", LUT))  LUT = lonlat2planar(LUT, proj.type=pL$aegis_proj4string_planar_km)
      
      if (! inherits(LOCS$timestamp, "POSIXct") )  LOCS$timestamp =  lubridate::date_decimal( LOCS$timestamp, tz=tz )
  
      if (!exists("dyear", LOCS) | (!exists("yr", LOCS)) ) {
        LOCS$tiyr = lubridate::decimal_date( LOCS$timestamp  ) 
        LOCS$yr = trunc( LOCS$tiyr )
        LOCS$dyear = LOCS$tiyr - LOCS$yr   
      } 

    
      LOCS$tplon = trunc(LOCS$plon / space_resolution + 1 ) * space_resolution
      LOCS$tplat = trunc(LOCS$plat / space_resolution + 1 ) * space_resolution
      LOCS$tdyear = trunc(LOCS$dyear / time_resolution + 1 ) * time_resolution

      nY = diff( range(LUT$yr) )
      nW = round(1/time_resolution)
      Tdims = c( nY, nW )
      Tres = c(1, time_resolution)
      Torigin = c( min(LUT$yr), 0) 

      plon_range = range( LUT$plon )
      plat_range = range( LUT$plat )
      plons = seq(plon_range[1], plon_range[2], by=space_resolution)
      plats = seq(plat_range[1], plat_range[2], by=space_resolution)
      Sdims = c(length(plons), length(plats))
      Sres = c(space_resolution, space_resolution)
      Sorigin = c( plon_range[1], plat_range[1] )

      ii = match( 
        paste(
          array_map( "xy->1", LOCS[, c("tplon","tplat")], dims=Sdims, res=Sres, origin=Sorigin ), 
          array_map( "ts->1", LOCS[, c("yr", "tdyear")], dims=Tdims, res=Tres, origin=Torigin ), 
          sep="_"
        ), 
        paste( 
          array_map( "xy->1", LUT[,c("plon","plat")], dims=Sdims, res=Sres, origin=Sorigin ), 
          array_map( "ts->1", LUT[,c("yr", "dyear")], dims=Tdims, res=Tres, origin=Torigin ), 
          sep="_" 
        ) 
      )

      LOCS$tplon = LOCS$tplat = LOCS$tdyear = NULL
      gc()

      for (vnm in variable_name) {
        if ( vnm %in% names(LUT )) {
          LOCS[[vnm]] = NA  
          LOCS[[vnm]] = LUT[ ii, vnm ]
        } 
      }



    }


    if ( project_class %in% c("core", "stmv", "hybrid") & output_format == "areal_units" )  {

      if (exists("lon", LUT)) {
        LUT = sf::st_as_sf( LUT, coords=c("lon", "lat") )
        st_crs(LUT) = st_crs( projection_proj4string("lonlat_wgs84") )
      } else if (exists("plon", LUT)) {
        LUT = sf::st_as_sf( LUT, coords=c("plon", "plat") )
        st_crs(LUT) = st_crs(pL$aegis_proj4string_planar_km) 
      } 
  
      Tdims = c( diff( range(LUT$yr) ), pL$nw )
      Tres = c(1, pL$tres)
      Torigin = c( min(LUT$yr), 0 ) 

      Sdims = c( pL$nplons, pL$nplats )
      Sres = c( pL$pres, pL$pres )
      Sorigin = pL$origin

      if (!exists("AUID", LOCS_AU))  {
        message ("AUID not found in LOCS_AU polygons, setting AUID as row number of polygons")
        LOCS_AU$AUID = as.character(1:nrow(LOCS_AU))
      }

      if ( ! "sf" %in% class(LOCS) ) {
        LOCS = st_as_sf( LOCS, coords=c("lon","lat"), crs=st_crs(projection_proj4string("lonlat_wgs84")) )
      }

      LOCS_AU = sf::st_transform( LOCS_AU, crs=st_crs(LUT) )

      if (! inherits(LOCS$timestamp, "POSIXct") )  LOCS$timestamp =  lubridate::date_decimal( LOCS$timestamp, tz=tz )
      if (! exists("yr", LOCS) )LOCS$yr = lubridate::year(LOCS$timestamp) 
      if (! exists("dyear", LOCS) ) LOCS$dyear = lubridate::decimal_date( LOCS$timestamp ) - LOCS$yr


      LU_map = paste( 
        st_points_in_polygons( pts=LUT, polys=LOCS_AU[, "AUID"], varname= "AUID" ), 
        array_map( "ts->1", st_drop_geometry( LUT) [,c("yr", "dyear")], dims=Tdims, res=Tres, origin=Torigin), 
        sep="_"
      )

      LOCS_map =  paste( 
        st_points_in_polygons( pts=LOCS, polys = LOCS_AU[, "AUID"], varname= "AUID" ),  
        array_map( "ts->1", st_drop_geometry(LOCS)[ , c("yr", "dyear") ], dims=Tdims, res=Tres, origin=Torigin), 
        sep="_"
      )


      LUT = setDT(st_drop_geometry(LUT))

      for (vnm in variable_name) {
        if ( vnm %in% names(LUT )) {  
          LUT_regridded = LUT[, setNames(.(FUNC( get(vnm), na.rm=TRUE) ), vnm), by=.(LU_map) ]
          LOCS[[ vnm ]] = LUT_regridded[ match( LOCS$AUID, LUT_regridded$LU_map ), vnm, with=FALSE ]
        }   
      }


    }



    if ( project_class %in% c("carstm" ) & output_format == "points" )  {
      # areal unit (LUT) to points (LOCS)
      nw = length(LUT$cyclic_name)
      ny = length(LUT$time_name)
      yr0 = min(as.numeric(LUT$time_name))

      if (is.null(LUT_AU)) LUT_AU = LUT$sppoly
      LUT_AU = st_make_valid(LUT_AU)

      LUT_AU = st_cast( LUT_AU, "MULTIPOLYGON")  # causes additional polys ..
      # LUT_AU = st_cast( LUT_AU, "POLYGON" )
      LUT_AU = st_make_valid(LUT_AU)
      LUT_AU = sf::st_transform( LUT_AU, crs=st_crs(pL$aegis_proj4string_planar_km) )
      LUT_AU$au_index = 1:nrow(LUT_AU)
      if (!exists("AUID", LUT_AU)) {
        message ("AUID not found in LOCS_AU polygons, setting AUID as row number of polygons")
        LUT_AU$AUID = as.character(LUT_AU$au_index)
      }
      bm = match( LUT_AU$AUID, LUT$space_name )  


      if (exists("lon", LOCS)) {
        LOCS = sf::st_as_sf( LOCS, coords=c("lon", "lat") )
        st_crs(LOCS) = st_crs( projection_proj4string("lonlat_wgs84") )
        LOCS = sf::st_transform( LOCS, crs=st_crs(pL$aegis_proj4string_planar_km) )
      } else if  (exists("plon", LOCS)) { 
        LOCS = sf::st_as_sf( LOCS, coords=c("lon", "lat") )
        st_crs(LOCS) = st_crs(pL$aegis_proj4string_planar_km)
      }

      LOCS$AUID = st_points_in_polygons( pts=LOCS, polys = LUT_AU[, "AUID"], varname= "AUID" )   

      if (! inherits(LOCS$timestamp, "POSIXct") )  LOCS$timestamp =  lubridate::date_decimal( LOCS$timestamp, tz=tz )
      if (! exists("yr", LOCS) ) LOCS$yr = lubridate::year(LOCS$timestamp) 
      if (! exists("dyear", LOCS) ) LOCS$dyear = lubridate::decimal_date( LOCS$timestamp ) - LOCS$yr
    

      ii = cbind( match( LOCS$AUID, LUT$space_name[bm] ) )
  
      jj = cbind( 
        ii, 
        array_map( "ts->year_index", LOCS[["yr"]], dims=c(ny ), res=c( 1  ), origin=c( yr0 ) ) 
      )


      LOCS_DF = LOCS
      if (inherits(LOCS_DF, "sf")) LOCS_DF = st_drop_geometry(LOCS_DF)

      ll = cbind( 
        ii, 
        array_map( "ts->2", LOCS_DF[, c("yr", "dyear")], dims=c(ny, nw), res=c( 1, 1/nw ), origin=c( yr0, 0) ) 
      )
      LOCS_DF = NULL

      for (g in 1:length(statvars)) {
        stat_var = statvars[g]
        for (h in 1:length(variable_name) ) {
          vnh = variable_name[[h]] 
          if (exists(vnh[[1]], LUT)) {
            vnm = paste( paste0(vnh,  collapse="_"),  stat_var, sep="_" )
            o = carstm_results_unpack( LUT, vnh ) 
            kk = length(dim(o))
            if (kk == 2) LOCS[[vnm]] = o[ cbind(ii, which(dimnames(o)$stat == stat_var) )  ] 
            if (kk == 3) LOCS[[vnm]] = o[ cbind(jj, which(dimnames(o)$stat == stat_var) )  ] 
            if (kk == 4) LOCS[[vnm]] = o[ cbind(ll, which(dimnames(o)$stat == stat_var) )  ] 
          }
        } 
      }
    
    }


    if ( project_class %in% c("carstm") & output_format == "areal_units" )  {
      # areal unit (LUT) to areal units (LUT_AU/LOCS) 
      # from source data: LUT = modelled predictions; LUT_AU are associated areal units linked by "AUID" 

      # format coordinate systems   
      if (is.null(LUT_AU)) LUT_AU = LUT$sppoly
      LUT_AU = st_make_valid(LUT_AU)
      LUT_AU = st_cast( LUT_AU, "MULTIPOLYGON")
      # LUT_AU = st_cast (LUT_AU, "POLYGON" )
      LUT_AU = st_make_valid(LUT_AU)
      LUT_AU = sf::st_transform( LUT_AU, crs=st_crs(pL$aegis_proj4string_planar_km) )

      nw = length(LUT$cyclic_name)
      ny = length(LUT$time_name)
      yr0 = min(as.numeric(LUT$time_name))

      if (!exists("AUID", LUT_AU)) LUT_AU$AUID = as.character(1:nrow(LUT_AU))
      bm = match( LUT_AU$AUID, LUT$space_name )  # should not be required but just in case things are subsetted

      # lut_uid is A LOCAL index of LUT_AU /LUT .. rasterize it
      LUT_AU$lut_uid = 1:nrow(LUT_AU)
    
      LUT_AU_raster = stars::st_rasterize( LUT_AU["lut_uid"], dx=space_resolution, dy=space_resolution )
      LUT_AU_pts = sf::st_as_sf( LUT_AU_raster, as_points=TRUE, na.rm=FALSE )
      st_crs(LUT_AU_pts) = st_crs( LUT_AU ) 

      # format output polygon structure 
      LOCS_AU = sf::st_transform( LOCS_AU, crs=st_crs(LUT_AU) )  # LOCS_AU .... must be sent ... <---------
      LOCS_AU = st_make_valid(LOCS_AU)
      LOCS_AU = st_buffer( LOCS_AU, 0 )
      LOCS_AU = st_cast( LOCS_AU, "MULTIPOLYGON")
      # LOCS_AU = st_cast( LOCS_AU, "POLYGON")
      if (!exists("AUID", LOCS_AU))  {
        message ("AUID not found in the LOCS_AU polygons, setting AUID as row number of polygons")
        LOCS_AU$AUID = as.character(1:nrow(LOCS_AU))
      }

      if (! inherits(LOCS, "sf") )   {
        if (exists("lon", LOCS )) {
          LOCS = st_as_sf( LOCS, coords=c("lon","lat"), crs=st_crs(projection_proj4string("lonlat_wgs84")) )
          LOCS = sf::st_transform( LOCS, crs=st_crs(pL$aegis_proj4string_planar_km) )
          if (!exists("AUID", LOCS )) LOCS[["AUID"]]  = st_points_in_polygons( pts=LOCS, polys=LOCS_AU[, "AUID"], varname= "AUID" ) 
        }
      }
      
      if (! inherits(LOCS$timestamp, "POSIXct") )  LOCS$timestamp =  lubridate::date_decimal( LOCS$timestamp, tz=tz )
      if (! exists("yr", LOCS) ) LOCS$yr = lubridate::year(LOCS$timestamp) 
      if (! exists("dyear", LOCS) ) LOCS$dyear = lubridate::decimal_date( LOCS$timestamp ) - LOCS$yr
  
      LOCS_DF = LOCS
      if (inherits(LOCS_DF, "sf")) LOCS_DF = st_drop_geometry(LOCS_DF)

      TIMESTAMP_index1 = array_map( "ts->year_index", LOCS_DF[, "yr"], dims=c(ny ), res=c( 1  ), origin=c( yr0 ) )
      TIMESTAMP_index2 = array_map( "ts->2", LOCS_DF[, c("yr", "dyear")], dims=c(ny, nw), res=c( 1, 1/nw ), origin=c( yr0, 0) )
      
      if (any( TIMESTAMP_index1 < 0)) {
        TIMESTAMP_index1[ which( TIMESTAMP_index1 <= 0 ) ] = NA
      }
      if (any( TIMESTAMP_index2 < 0)) {
        warning ("cyclic time index is negative: lookup table does not contain data to lookup ... \n 
          lookup data needs to be expanded or time periods reduced in model")
        TIMESTAMP_index2[ which( TIMESTAMP_index2 <= 0 ) ] = NA
      }

      LOCS_DF = NULL

      # id membership of LUT raster in LOCS_AU
      LUT_AU_pts_LOCS_AU_AUID = st_points_in_polygons( pts=LUT_AU_pts, polys=LOCS_AU[,"AUID"], varname="AUID" ) 

      pts_AU = match(LUT_AU_pts$lut_uid, LUT_AU$lut_uid[match(LUT$space_name, LUT_AU$AUID)] ) ## (layer==lut_uid) -- to -- LUT
      aggFUN = function( LUV ) {
          tapply(X=LUV, INDEX=LUT_AU_pts_LOCS_AU_AUID, FUN=FUNC, na.rm=TRUE) 
      }

      for (g in 1:length(statvars)) {
        stat_var = statvars[g]
        for (h in 1:length(variable_name) ) {
          vnh = variable_name[[h]] 
          if (exists(vnh[[1]], LUT)) {
            vnm = paste( paste0(vnh,  collapse="_"),  stat_var, sep="_" )
            o = carstm_results_unpack( LUT, vnh ) 
            kk = length(dim(o))
            if (kk == 2) {
              ov = o[ pts_AU,  which(dimnames(o)$stat == stat_var)  ] 
              LUT_as_LOCS_AU = aggFUN( ov)
              space_index = match( as.character(LOCS$AUID), as.character(dimnames(LUT_as_LOCS_AU )[[1]] )  )   # AUID of LOCS_AU (from LUT_AU_pts_AUID)
              if (all(is.na(space_index))) {
                LOCS[,vnm] = LUT_as_LOCS_AU  # fiddling required when only 1 AU
              } else {
                LOCS[,vnm] = LUT_as_LOCS_AU[  space_index  ]
              }
            }
            if (kk == 3) {
              ov = o[ pts_AU,, which(dimnames(o)$stat == stat_var)  ] 
              LUT_as_LOCS_AU = apply( ov, MARGIN=c(2), FUN=aggFUN )
              space_index = match( as.character(LOCS$AUID), as.character(dimnames(LUT_as_LOCS_AU )[[1]] )  )   # AUID of LOCS_AU (from LUT_AU_pts_AUID)
              if (all(is.na(space_index))) {
                LOCS[,vnm] = LUT_as_LOCS_AU[ TIMESTAMP_index1  ]  # fiddling required when only 1 AU
              } else {
                LOCS[,vnm] = LUT_as_LOCS_AU[ cbind( space_index, TIMESTAMP_index1 ) ]
              }
            }
            if (kk == 4) {
              ov = o[ pts_AU,,, which(dimnames(o)$stat == stat_var)  ] 
              LUT_as_LOCS_AU = apply( ov, MARGIN=c(2,3), FUN=aggFUN )
              space_index = match( as.character(LOCS$AUID), as.character(dimnames(LUT_as_LOCS_AU )[[1]] )  )   # AUID of LOCS_AU (from LUT_AU_pts_AUID)
              if (all(is.na(space_index))) {
                LOCS[,vnm] = LUT_as_LOCS_AU[ TIMESTAMP_index2  ]  # fiddling required when only 1 AU
              } else {
                LOCS[,vnm] = LUT_as_LOCS_AU[ cbind( space_index, TIMESTAMP_index2 ) ]
              }
            }
          }
        } 
      }
    }
  } # end space-time-cyclic


  if (returntype =="sf" ) {
    if (! inherits(LOCS, "sf"))   LOCS = st_as_sf( LOCS, coords=c("lon","lat"), crs=st_crs(projection_proj4string("lonlat_wgs84")) )
  }
  if (returntype =="data.frame" ) {
    if  (! inherits(LOCS, "data.frame"))  LOCS = as.data.frame(LOCS)
  }

  if (returntype =="data.table" ) {
    if (! inherits(LOCS, "data.table"))   setDT(LOCS)
  }

  if (returntype =="vector" ) {
    if (length(variable_name) == 1) {
      if (project_class=="carstm") {
        variable_name = paste( paste0(variable_name[[1]],  collapse="_"),  statvars[1], sep="_" )
      }          
      if (exists( variable_name, LOCS)) {
        LOCS = LOCS[[variable_name]] 
      } else {
        LOCS = rep(NA, nrow(LOCS))
      }
    }  
  }

  return(LOCS)

}
jae0/aegis documentation built on Jan. 14, 2025, 7:46 p.m.