Extent Indicators, Zooming and Clearer Maps with grmap

[Jared Greathouse]

I'm doing a synthetic control analysis of COVID-19 vaccine mandates in the United States, and I want to use grmap to show how treatment and donor pool counties are imbalanced on key background covariates (e.g., urbanicity, election results, comorbidities). I want to map a given region and then map a subset of said region alongside it. I found two similar posts on Statalist- one here, and another here- yet I cannot use these as solutions because little to no code was given.

My analysis divides the U.S. into regions. Here's my code so you can follow me.

cd "D:\Test" // Or whatever directory you like


cap conf f us_final.dta

if _rc { // if file doesn't exist, as is likely

loc file cb_2018_us_county_500k.zip


copy "https://www2.census.gov/geo/tiger/GENZ2018/shp/`file'" ///
 "`file'", replace
 

unzipfile `file', replace


* Makes Stata format shapefile
spshape2dta cb_2018_us_county_500k.dbf, replace saving(us_final)

conf f us_file

foreach name in cb {

local files : dir "`c(pwd)'" files "*`name'*"

foreach f of loc files { // erases old files
    
erase  `f'
}
}
}

u us_final, clear

cls // clears the screen

rename (STATEFP COUNTYFP) (sid cid) // state id, county id

destring *id, replace

/*

Making the regions

*/

g region = .

replace region = 1 if inlist(sid,9,10,11,23,24,25,33,34,36,42,44,50,51) // Northeast Region

replace region = 2 if inlist(sid,1,5,12,13,21,22,28,37,45,47,48,54) // South


replace region = 3 if inlist(sid,6,41,53) // Pacific West

keep if inrange(region,1,3)
as !mi(region)


grmap if region ==1

Okay, so we see the Northeastern U.S. I like the regional map, but my treatment area for region 1 is New York City, or Spatial IDs 2835, 170, 990, 2218 and 166.

To iterate, I just want to take these spatial IDs and put them in an extent indicator/box to the side and zoom in on it, so I and my readers can have better vision of my treatment area compared to its donor pool. To give a visual example of what exactly I mean, please see this post or this one on GIS stack exchange/overflow. Is this possible in Stata using grmap, or must I resort to a user-written command?

[Anne Denis]

Hello,

out of curiosity did you find a way to do this? I need to do the exact same thing, and I know how to do this in other languages such as python but have to use STATA, and I don't even know it is possible to do so? Currently I am using spmap, and have yet to find a solution

[Ali Atia]

Some basic technique (this code should replace the last line of code in #1):

grmap if inlist(_ID,2835,170,990,2218,166),name(zoomedin) fysize(25) fxsize(25) plotregion(lcolor(black))
sum _CX,detail
local xmax = r(max)
preserve
keep if inlist(_ID,2835)
rename (_CX _CY) (_X1 _Y1)
gen _X2 = `xmax'
gen _Y2 = 41
save arrow.dta,replace
restore
grmap if region ==1, arrow(data(arrow)) name(full) point(data(us_final) select(keep if inlist(_ID,2835)) xcoord(_CX) ycoord(_CY) fcolor(red) ocolor(red))
gr combine full zoomedin, commonscheme scheme(s1mono)

[Anne Denis]

Hi Ali,

Thank you very much for your response, this is exactly what I am looking for. However, I am trying to find more information on grmap--but stata says it is unrecognized, and I am unable to install it with ssc. the only iformation I have found on it is this: https://www.stata.com/manuals/spgrmap.pdf , but unfortunately it does not provide any information on syntax. I feel like I am missing something obvious here but I'm not sure what (I am new to mapping with stata)! In any case, thank you so much for your response it is much appreciated!!

[Ali Atia]

You need to activate grmap before you can use it:

 grmap, activate

[Jared Greathouse]

Ali Atia This is so cool, thanks a lot, I didn't even know if there was a way to do this.