search

aknandi / disaggregation

R package containing methods for Bayesian disaggregation modelling
Other
21 stars 6 forks source link

predict_iid = TRUE in predict_uncertainty doesn't work in my example. #26

Closed timcdlucas closed 5 years ago

timcdlucas commented 5 years ago 

library(disaggregation)
library(raster)
library(dplyr)

source("setUserInfo.R")

# define paths

PR_path <- Z('GBD2017/Processing/Stages/04b_PR_DB_Import_Export/Verified_Outputs/2018_02_15/pfpr_dump.csv')
API_path <- Z('GBD2017/Processing/Stages/04c_API_Data_Export/Checkpoint_Outputs/subnational.csv')
pop_path <- Z('GBD2017/Processing/Stages/03_Muster_Population_Figures/Verified_Outputs/Output_Pop_At_Risk_Pf_5K/ihme_corrected_frankenpop_All_Ages_3_2015_at_risk_pf.tif')
shapefile_path <- Z('master_geometries/Admin_Units/Global/GBD/GBD2017_MAP/GBD2017_MAP_MG_5K/MG_5K_GBD2017_ADMIN3.shp')

cov_raster_paths <- c(
  Z('mastergrids/MODIS_Global/MOD11A2_v6_LST/LST_Day/5km/Synoptic/LST_Day_v6.Synoptic.Overall.mean.5km.mean.tif'),
  Z('mastergrids/MODIS_Global/MCD43D6_v6_BRDF_Reflectance/EVI_v6/5km/Synoptic/EVI_v6.Synoptic.Overall.mean.5km.mean.tif'),
  #Z('mastergrids/Other_Global_Covariates/TemperatureSuitability/TSI_Pf_Dynamic/5km/Synoptic/TSI-Martens2-Pf.Synoptic.Overall.Mean.5km.Data.tif'),
  #Z('GBD2017/Processing/Static_Covariates/MAP/other_rasters/accessibility/accessibility.5k.MEAN.tif'),
  Z('mastergrids/Other_Global_Covariates/Elevation/SRTM-Elevation/5km/Synoptic/SRTM_elevation.Synoptic.Overall.Data.5km.mean.tif'),
  Z('mastergrids/MODIS_Global/MOD11A2_v6_LST/LST_Day/5km/Synoptic/LST_Day_v6.Synoptic.Overall.SD.5km.mean.tif'),
  #Z('mastergrids/MODIS_Global/MCD43B4_BRDF_Reflectance/TCB/5km/Synoptic/TCB.Synoptic.Overall.mean.5km.mean.tif'),
  #Z('mastergrids/Other_Global_Covariates/NightTimeLights/VIIRS_DNB_Composites/5km/Annual/VIIRS-SLC.2016.Annual.mean.5km.median.tif'),
  #Z('mastergrids/Other_Global_Covariates/UrbanAreas/Global_Urban_Footprint/From_86m/5km/Global_Urban_Footprint_5km_PropUrban.tif'),
  Z('mastergrids/MODIS_Global/MCD43D6_v6_BRDF_Reflectance/TCW_v6/5km/Synoptic/TCW_v6.Synoptic.Overall.mean.5km.mean.tif')
)

shps <- shapefile(shapefile_path)

covs <- lapply(cov_raster_paths, raster)
cov_stack <- stack(covs)

cov_stack <- crop(cov_stack, extent(c(40, 52, -26, -10)))

api <- read.csv(API_path)
api <-
  api %>%
    filter(iso3 == 'MDG', year == 2013)

shps <- shps[shps$area_id %in% as.character(api$shapefile_id), ]

shps$inc <- api$api_mean_pf[match(shps$area_id, as.character(api$shapefile_id))]

pop_ras <- raster(pop_path)
pop_ras <- crop(pop_ras, extent(c(40, 52, -26, -10)))

cov_stack <- crop(cov_stack, pop_ras)

dis_data <- prepare_data(shps, cov_stack, 
                         aggregation_raster = pop_ras, 
                         id_var ='area_id',
                         response_var = 'inc',
                         ncores = 8)

dis_data$covariate_data[is.na(dis_data$covariate_data)] <- 0

dis_data$aggregation_pixels[is.na(dis_data$aggregation_pixels)] <- 0

m <- fit_model(dis_data, family = 'gaussian', link = 'identity', field = FALSE)

mpred <- predict_model(m)
munc <- predict_uncertainty(m, N = 40, predict_iid = F)
munc2 <- predict_uncertainty(m, N = 40, predict_iid = T)
 Error in p$rasterize(nrow(r), ncol(r), as.vector(extent(r)), values, background) :

Not very minimal. I'm going to look at this. But thought I'd document it here in case I get pulled away.