maps for Year 1 report

seattleflu · Jan 12, 2020 · 8eb2c9c · 8eb2c9c
1 parent 35b9714
commit 8eb2c9c
Show file tree

Hide file tree

Showing 3 changed files with 46 additions and 19 deletions.
diff --git a/exploration_scripts/model_maps/Year1_catchment.R b/exploration_scripts/model_maps/Year1_catchment.R
@@ -18,8 +18,8 @@ g_legend <- function(a.gplot){
 } 
 
 # setting
-SOURCE='seattle_geojson'
-GEO = 'residence_neighborhood_district_name'
+SOURCE='sfs_domain_geojson'
+GEO = 'residence_regional_name'
 
 
 # build model
@@ -35,11 +35,13 @@ db <- expandDB(selectFromDB(  queryIn, source='production', na.rm=TRUE ), shp=sh
 
 ## set up map environment
 plotDat <- right_join(db$observedData,shp, by=GEO)
-plotDat <- plotDat %>% group_by(site_type) %>% mutate(fraction = n/max(n))
+plotDat <- plotDat %>% group_by(site_type) %>% mutate(fraction = n/max(n)) %>%
+  filter(site_type %in% c('childcare','clinic','collegeCampus','homelessShelter',
+                                        'port','retrospective','self-test','workplace'))
 
 bbox<-sf::st_bbox(shp$geometry)
 
-mapSettings <- ggplot() + #xlim(c(min(122.5, -bbox[1]),max(121.7,-bbox[3]))) + ylim(c(max(47.17,bbox[2]),min(47.76,bbox[4]))) +
+mapSettings <- ggplot() + xlim(-122.45,-122.05) + ylim(47.45,47.8) +
   theme_bw() +
   theme(axis.text=element_blank(),axis.ticks=element_blank(),panel.grid.major=element_line(colour="transparent"), panel.border = element_blank()) +
   geom_sf(data=shp,size=0.1,aes(fill=NaN)) 
@@ -48,18 +50,24 @@ mapSettings <- ggplot() + #xlim(c(min(122.5, -bbox[1]),max(121.7,-bbox[3]))) + y
 # colorLimits<-c(min(plotDat$n,na.rm=TRUE),max(plotDat$n,na.rm=TRUE))
 # colorBreaks<-unique(c(min(colorLimits),round(1e0*seq(sqrt(min(colorLimits)),sqrt(max(colorLimits)), length.out = 5)^2)/1e0, max(colorLimits)))
 
+plotDat$site_type[plotDat$site_type=='clinic']<-"childrensOutpatient"
+plotDat$site_type[plotDat$site_type=='port']<-"Sea-Tac"
+plotDat$site_type[plotDat$site_type=='retrospective']<-"hospitalRetrospective"
+
+
 p <- mapSettings + geom_sf(data=plotDat ,size=0, aes(fill=fraction, group='site_type'))  +
-  guides(fill=guide_legend(title="Year 1\n relative recruitment")) +
+  guides(fill=guide_legend(title="  relative\nrecruitment")) +
   viridis::scale_fill_viridis(na.value="transparent") + #,breaks=colorBreaks,limits=colorLimits)+ 
-  facet_wrap('site_type')
+  facet_wrap('site_type', nrow=2)
 p
 
-png(filename=paste('Year1_catchment_map_by_site_type.png',sep='-'),res=600, units='in', width=10,height=6)
+png(filename=paste('Year1_catchment_map_by_site_type.png',sep='-'),res=600, units='in', width=6.67,height=4)
 p
 dev.off()
 
 
-plotDat2 <- plotDat %>% group_by(!!rlang::parse_expr(GEO)) %>% summarize(n=sum(n)) %>% mutate(fraction = n/max(n)) %>% right_join(shp,by=GEO)
+plotDat2 <- plotDat %>% group_by(!!rlang::parse_expr(GEO)) %>% summarize(n=sum(n)) %>% mutate(fraction = n/max(n)) %>% right_join(shp,by=GEO) 
+
 p2 <- mapSettings + geom_sf(data=plotDat2 ,size=0, aes(fill=n))  +
   guides(fill=guide_legend(title="Year 1 count")) +
   viridis::scale_fill_viridis(na.value="transparent",breaks=seq(350,1600,by=200))

diff --git a/exploration_scripts/model_maps/Year1_catchment_map_by_site_type.png b/exploration_scripts/model_maps/Year1_catchment_map_by_site_type.png
diff --git a/exploration_scripts/model_maps/Year1_fourMaps.R b/exploration_scripts/model_maps/Year1_fourMaps.R
@@ -23,25 +23,32 @@ fluPathogens <- c('Flu_A_H1','Flu_A_H3','Flu_A_pan','Flu_B_pan','Flu_C_pan')
 SOURCE='sfs_domain_geojson'
 GEO = 'residence_regional_name'
 
+
+# SOURCE='seattle_geojson'
+# GEO = 'residence_neighborhood_district_name'
+
+factors   <- c('site_type')#,'age_range_fine_upper')
+shp <- masterSpatialDB(shape_level = gsub('residence_','',GEO), source = SOURCE)
+
 # SOURCE='seattle_geojson'
 # GEO = 'residence_neighborhood_district_name'
 
-PATHOGEN='Flu_A_H1'
+# PATHOGEN='Flu_A_H1'
 # PATHOGEN='Flu_A_H3'
 # PATHOGEN <- fluPathogens
 # PATHOGEN='RSVA'
 # PATHOGEN='RSVB'
 # PATHOGEN=c('RSVA','RSVB')
+PATHOGEN='all'
+# PATHOGEN='CoV'
 
 # db <- selectFromDB(queryIn= list(SELECT  =c("*")), source = SRC)
 # pathogens <- db$observedData %>% group_by(pathogen) %>% summarize(n = n()) %>% arrange(desc(n))
 # PATHOGEN = setdiff(pathogens$pathogen,c(fluPathogens,'not_yet_tested','measles','Measles'))
 
 
-factors   <- c('site_type','sex','flu_shot')#,'age_range_fine_upper')
 
 # build model
-shp <- masterSpatialDB(shape_level = gsub('residence_','',GEO), source = SOURCE)
 queryIn <- list(
   SELECT   =list(COLUMN=c('pathogen', factors, GEO,'encountered_week')),
   WHERE    =list(COLUMN='pathogen', IN=PATHOGEN),
@@ -50,11 +57,11 @@ queryIn <- list(
 )
 
 db <- expandDB(selectFromDB(  queryIn, source='production', na.rm=TRUE ), shp=shp)
-
 db <- appendCatchmentModel(db,shp=shp, source='production', na.rm=TRUE )
-modelDefinition <- latentFieldModel(db=db, shp=shp)
 
+modelDefinition <- latentFieldModel(db=db, shp=shp)
 model <- modelTrainR(modelDefinition)
+summary(model$inla)
 
 ## line plot
 ggplot(model$latentField) + 
@@ -69,28 +76,40 @@ ggplot(model$latentField) +
 plotDat <- right_join(model$latentField %>% group_by_(.dots =c(GEO,'encountered_week')) %>% summarise(modeled_intensity_median = median(modeled_intensity_median)),shp, by=GEO)
 bbox<-sf::st_bbox(shp$geometry)
 
-mapSettings <- ggplot() + xlim(-122.45,-122.05) + ylim(47.3,47.75) +
+# #h1
+# plotDat$modeled_intensity_median <- pmin(0.10,plotDat$modeled_intensity_median)/0.10
+# 
+# #h3
+# plotDat$modeled_intensity_median <- pmin(0.12,plotDat$modeled_intensity_median)/0.12
+
+#all
+plotDat$modeled_intensity_median <- pmin(1,plotDat$modeled_intensity_median)/1
+
+mapSettings <- ggplot() + xlim(-122.45,-122.05) + ylim(47.45,47.8) +
   theme_bw() +
   theme(axis.text=element_blank(),axis.ticks=element_blank(),panel.grid.major=element_line(colour="transparent"), panel.border = element_blank()) +
   geom_sf(data=shp,size=0.1,aes(fill=NaN))
 
-colorLimits<-c(min(plotDat$modeled_intensity_median,na.rm=TRUE),max(plotDat$modeled_intensity_median,na.rm=TRUE))
-colorBreaks<-unique(c(min(colorLimits),round(1e4*seq(sqrt(min(colorLimits)),sqrt(max(colorLimits)), length.out = 5)^2)/1e4, max(colorLimits)))
+colorLimits=c(0,1)
+colorBreaks=seq(0,5,by=1)/5
+# colorLimits<-c(min(plotDat$modeled_intensity_median,na.rm=TRUE),max(plotDat$modeled_intensity_median,na.rm=TRUE))
+# colorBreaks<-unique(c(min(colorLimits),round(1e4*seq(sqrt(min(colorLimits)),sqrt(max(colorLimits)), length.out = 5)^2)/1e4, max(colorLimits)))
 
 ## choose weeks to show
-week <- c('2019-W01','2019-W05','2019-W09', '2019-W13','2019-W17')
+week <- c('2019-W02','2019-W05','2019-W08', '2019-W11','2019-W14')
 p <- list()
 for(k in 1:length(week)){
   p[[k]] <- mapSettings + geom_sf(data=plotDat %>% filter(encountered_week == week[k]),size=0, aes(fill=modeled_intensity_median))  +
     guides(fill=guide_legend(title="modeled_intensity")) +
-    viridis::scale_fill_viridis(na.value="transparent",trans = "sqrt",breaks=colorBreaks,limits=colorLimits)+
+    viridis::scale_fill_viridis(na.value="transparent",breaks=colorBreaks,limits=colorLimits)+
     ggtitle(week[k])
   if(k<length(week)){
     p[[k]] <- p[[k]] + theme(legend.position='none')
   }
 }
 legend <- g_legend(p[[k]]) 
+p[[k]] <- p[[k]] + theme(legend.position='none')
 
 png(filename=paste('Year1',paste(PATHOGEN,collapse='-',sep=''),'fourMaps.png',sep='-'),res=600, units='in', width=6.5,height=4)
-grid.arrange(p[[1]],p[[2]],p[[3]],p[[4]], legend,nrow=1)
+grid.arrange(p[[1]],p[[2]],p[[3]],p[[4]], p[[5]], legend,nrow=1)
 dev.off()