I am working on mapping an outcome variable (% of mothers receiving postnatal care) at the county level in Kenya using srmap and shp2dta both from SSC in Stata 17 on a Mac. I was able to create the map but am struggling with two issues: (1) getting labels with the county name onto the map and (2) making the counties with *low* postnatal care coverage darker. It looks like all the sequential color schemes in Stata run darker with higher values, which usually works but I want to highlight the low coverage areas. I can’t use a divergent scheme because I need all counties to have a color. Looks like a reverse option is not allowed, any ideas on this? I want to present % who received pnc, not % who didn’t so I’m not able to flip the variable as a workaround.
The more important issue though is getting the county names onto the map. Looks like maybe I need to create another dataset with just the county names and x,y coordinates? Getting confused though, would appreciate any help.
Below is the code that I have so far. Note: I’m sorry for the inconvenience, but for the below code to run, the shapefiles would have to be downloaded from the internet to your computer and filepath below updated. The shapefiles are in a zipfile (ken_adm_iebc_20191031_SHP.zip) (second link down) at this website: https://data.humdata.org/dataset/cod-ab-ken. I tried to write code to read these shapefiles directly into STATA and unzip for a reproducible example, but I’m not sure it’s possible since there is only a download option and not a unique URL...
Below is the code I have so far. Thank you very much in advance for any help!
The more important issue though is getting the county names onto the map. Looks like maybe I need to create another dataset with just the county names and x,y coordinates? Getting confused though, would appreciate any help.
Below is the code that I have so far. Note: I’m sorry for the inconvenience, but for the below code to run, the shapefiles would have to be downloaded from the internet to your computer and filepath below updated. The shapefiles are in a zipfile (ken_adm_iebc_20191031_SHP.zip) (second link down) at this website: https://data.humdata.org/dataset/cod-ab-ken. I tried to write code to read these shapefiles directly into STATA and unzip for a reproducible example, but I’m not sure it’s possible since there is only a download option and not a unique URL...
Below is the code I have so far. Thank you very much in advance for any help!
Code:
*Generate example dataset
clear
set obs 20
gen id = _n
set seed 1234
gen byte pnc = (runiform() < 0.75)
gen county = ""
replace county = "Kilifi" in 1/2
replace county = "Meru" in 3/5
replace county = "Nyeri" in 5/10
replace county = "Kajiado" in 11/15
replace county = "Siaya" in 16/20
lab define pnc 0 "No" 1 "Yes"
lab values pnc pnc
*Create dataset with county-level average PNC
preserve
sort county
collapse (mean) pnc, by(county)
rename county ADM1_EN
gen pnc_county = round(pnc * 100, 0.01)
save pnc_county_example_formap, replace
restore
*Create county-level map of pnc
shp2dta using "$controls/Maps/HUM_Kenya_adm_20Feb2020/ken_admbnda_adm1_iebc_20191031", ///
database(kenyadb) coordinates(kenyacoord) genid(id) replace
use kenyadb, clear
merge m:1 ADM1_EN using pnc_county_example_formap
spmap pnc_county using kenyacoord, id(id) fcolor(Reds2) ///
clnumber(3) ocolor(Black) osize(medium) ///
title ("Percent receiving postnatal care by study county", size(small)) ///
legtitle("% receiving maternal PNC ") ///
legorder(hilo) ///
legstyle(2) legend(region(lcolor(dknavy))) ///
plotregion(icolor(gs15)) graphregion(icolor(gs15))
Comment