Dear Statalisters,
I'm curious if anyone has any good ideas about how to smooth map
polygons in Stata. Most of the shape files out there on the web seem
to be very high-resolution, perhaps suitable for navigation but less
suitable for conveying rough distributions in a graphic (preferably
with a small file size). My approach in the past has been to drop
smaller polygons, pick a random subset of points in all remaining
polygons, and then to march across observations to put back in points
that seem to measure large changes in position. This approach is
illustrated below.
clear all
cap copy http://pped.org/stata/uscoord.dta /uscoord.dta
use /uscoord.dta, clear
*parameters tol and mod determine coarseness of results
local tol 1
local mod 200
g c=sum(mi(_X))
g n=_n
egen count=count(_X), by(_I c)
g nc=-count
*keep only the important polygons
bys _ID (nc n): g p=sum(mi(_X))
keep if (p<2) | (p<3 & _ID==24) | (p<5 & _ID==13)
*get a fraction of points on each polygon
by _ID nc: g y=_Y if mod(_n,int(_N/`mod'))==2 | _ID==1
g x=_X if y<.
sort _ID n
*add back in far-distant points on each polygon
local r=6371/(2.54*5.280*.12)
loc p1 "sin(_X*_pi/180)*sin(_X[_n-1]*_pi/180)+cos(_X*_pi/180)"
loc p2 "*cos(_X[_n-1]*_pi/180)*cos((_Y[_n-1]-_Y)*_pi/180)"
g d=acos(`p1'`p2')*`r'
replace y=_Y if (d>`tol' & d<.)|(d[_n+1]>`tol'& d[_n+1]<.)
replace x=_X if (d>`tol' & d<.)|(d[_n+1]>`tol'& d[_n+1]<.)
drop if mi(y) & !mi(_Y)
*drop all but 50 states and DC
drop if inlist(_ID, 7, 40, 48, 54, 55)
* doctor loc and size of AK, HI, DC
replace y=y/2+10 if _ID==56
replace x=x/3-85 if _ID==56
replace x=x*1.5+102 if _ID==13
replace y=y+10 if _ID==13
replace y=_Y*10-355 if _ID==1
replace x=_X*10+700 if _ID==1
*drop extraneous obs
bys _ID p (n): drop if mi(y)&mi(y[_n-1])&_n>1
*make an obs to connect first and last points
bys _ID p (n): g expand=1+(_n==_N)
expand expand
bys _ID p (n): replace n=n[2]-1/2 if _n==_N
sort _ID n
loc o1 "name(usa, replace) xla(none) yla(none, nogrid) "
loc o2 "xsize(4) ysize(2) graphr(fc(white))xsc(noline) "
loc o3 `"ysc(noline) xti("") yti("") leg(off) "'
sc y x, c(l) ms(none) cmissing(n) `o1' `o2' `o3'
su x, meanonly
di %3.2f r(N)/521238*100 "% of orig data"
replace _X=x
replace _Y=y
keep _ID _X _Y
save /us, replace
use http://pped.org/stata/spop90.dta, clear
cap spmap pop using /us, id(id)
if _rc {
ssc inst spmap
spmap pop using /us, id(id)
}
If anyone has any better ideas, please let me know.
Thanks,
Austin
*
* For searches and help try:
* http://www.stata.com/help.cgi?search
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
