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.

Code:
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?