How to Create map of several U.S. with 3 circles and also use proc ginside to show points

[Juno77]

I need to create a map with 3 circles with different distances from the center point, but also to add markers for cities within each of the 3 circles. The code for the 3 circles covering a few states is below - this is thanks to an example I found online - robslink.com. The problem I have is incorporating the dataset with cities and using proc ginside to determine if the points are within a circle. Sincere Thank you to anyone who can help!

goptions reset=all cback=white border htitle=12pt htext=10pt xpixels=1000 ypixels=1000*/;
data dodge1;
input srvarea long lat city $25.;
x=atan(1)/45 * long;
y=atan(1)/45 * lat;
datalines;
100 100.0171 37.7528 Dodge City
;
run;
data dodge2;
input srvarea long lat city $25.;
x=atan(1)/45 * long;
y=atan(1)/45 * lat;

datalines;
200 100.0171 37.7528 Dodge City
;
run;
data dodge3;
input srvarea long lat city $25.;
x=atan(1)/45 * long;
y=atan(1)/45 * lat;

datalines;
300 100.0171 37.7528 Dodge City
;
run;
/********************************************/
data maplabel;
length function $ 8;
retain flag 0 xsys ysys '2' hsys '3' when 'a' style "'Albany AMT'";
set maps.uscenter(drop=long lat);

where (fipstate(state) in ('CO' 'NE' 'MO' 'NM' 'TX' 'KS' 'IA' 'OK'));
function='label'; text=fipstate(state); size=2.5; position='10';
if ocean='Y' then
do;
position='6'; output;
function='move';
flag=1;
end;

else if flag=1 then
do;
function='draw'; size=.5;
flag=0;
end;
output;
run;

/********************CIRCLE 1 100 MILES ***************/
data anno1;
set dodge1;
retain xsys ysys '2' flag 1 when 'a';
length text $25 color function $ 8 style $ 25;
drop xold yold;

xold=x;
yold=y;

d2r=3.1415926/180;
r=3958.739565;

xcen=long;
ycen=lat;

do degree=0 to 360 by 5;

if degree=0 then do;

function='poly';
line=1;
pattern v=e;
color='aF8A10266';
end;

else do;
function='polycont';

color='black';
end;
y=arsin(cos(degree*d2r)*sin(srvarea/R)*cos(ycen*d2r)+ cos(srvarea/R)*sin(ycen*d2r))/d2r;
x=xcen+arsin(sin(degree*d2r)*sin(srvarea/R)/cos(y*d2r))/d2r;

x=atan(1)/45*x;
y=atan(1)/45*y;
output;
end;
x=xold;
y=yold;

function='label';
style='special';
text='J';
position='5';
color='black';
output;

style="'Albany AMT'";
text=city;
position='2';
color='black';
output;
text=trim(left(srvarea))||' Miles';
size=1;
position='8';
output;

run;
//********************CIRCLE 2 - 200 MILES **************//
data anno2;
set dodge2;
retain xsys ysys '2' flag 1 when 'a';
length text $25 color function $ 8 style $ 25;
drop xold yold;

xold=x;
yold=y;
d2r=3.1415926/180;
r=3958.739565;
xcen=long;
ycen=lat;
do degree=0 to 360 by 5;
if degree=0 then do;
function='poly';
line=1;
pattern v=e;
color='aF8A10266';
end;
else do;
function='polycont';
color='blue';
end;
y=arsin(cos(degree*d2r)*sin(srvarea/R)*cos(ycen*d2r)+ cos(srvarea/R)*sin(ycen*d2r))/d2r;
x=xcen+arsin(sin(degree*d2r)*sin(srvarea/R)/cos(y*d2r))/d2r;
/* Convert degrees to radians */
x=atan(1)/45*x;
y=atan(1)/45*y;
output;
end;
x=xold;
y=yold;
function='label';
style='special';
text='J';
position='5';
color='blue';
output;
style="'Albany AMT'";
text=city;
position='2';
color='blue';
output;
text=trim(left(srvarea))||' Miles';
size=1;
position='8';
output;
run;
/***************** CIRCLE 3 - 300 MILES **************************/
data anno3;
set dodge3;
retain xsys ysys '2' flag 1 when 'a';
length text $25 color function $ 8 style $ 25;
drop xold yold;

xold=x;
yold=y;
d2r=3.1415926/180;
r=3958.739565;
xcen=long;
ycen=lat;
do degree=0 to 360 by 5;
if degree=0 then do;
function='poly';
line=1;
/*style='solid'; */
pattern v=e;
color='aF8A10266';
end;
else do;
function='polycont';
color='red';
end;
y=arsin(cos(degree*d2r)*sin(srvarea/R)*cos(ycen*d2r)+ cos(srvarea/R)*sin(ycen*d2r))/d2r;
x=xcen+arsin(sin(degree*d2r)*sin(srvarea/R)/cos(y*d2r))/d2r;
x=atan(1)/45*x;
y=atan(1)/45*y;
output;
end;
x=xold;
y=yold;
function='label';
style='special';
text='J';
position='5';
color='red';
output;
style="'Albany AMT'";
text=city;
position='2';
color='red';
output;

text=trim(left(srvarea))||' Miles';
size=1;
position='8';
output;
run;
/*** OTHER ANNOTATIONS ******/
data anno4;
function='label'; text=fipstate(state); size=2.5; position='10';
run;
data anno5;
set maplabel label1 maplabel;
run;
pattern1 r=52 v=msolid c=grayEE;
title1 '';
proc gmap data=ks map=ks anno=anno5;
id state;
choro state / coutline=black nolegend;
run;
quit;

CODE I TRIED FOR PROC GINSIDE TO INCLUDE POINTS INSIDE MAP CIRCLES - HAS PROBLEMS

data circle_map; set anno1;
inside_id=67801;
run;
/* proc ginside is a new v9.2 proc */
proc ginside data=kspk map=circle_map out=kspk;
id inside_id;
run;

/* Now, for dots that are 'inside' the circle, set their color to something special */
data kspk; set kspk
if (inside_id ne .) then color="red";
run;

/* Get the NC map */
proc sql;
create table mymap as
select state, segment, x, y
from maps.states
where
quit;
run;

data combined;
set mymap customers circle;
run;
proc gproject data=combined out=combined dupok ;
id state;
run;
data mymap kspk circle;
set combined;
if anno_flag=1 then output kspk;
if anno_flag=2 then output circle;
else output mymap;
run;

data anno_pk;
retain xsys ysys "2" hsys "3" when "a";
set kspk;
function="symbol"; style="marker"; text="V"; color="blue"; size=1.75;
output;
run;
/* Combine the annotate data set with the map data set */
data combo1;
set anno1 anno_pk maps.states(where=(state in (8 19 20 29 31 35 40 48)));
run;

[GraphGuy]

This example isn't exactly what you're wanting to do, but it should have several bits of code you can re-purpose ...

http://robslink.com/SAS/democd29/customer_info.htm

[Juno77]

Thank you for your response. I'm still not clear as to what the statement "id inside_id" does. I'm using cities, not zipcodes, do I have to use zipcodes to use proc ginside? How do I incorporate the cities annotation into all the other annotations. I tried it already by adding to the final anno ___ for proc gmap and it's not working.

Thanks again

[GraphGuy · Accepted solution]

Read through the explanation in my book (Example 9, starting on p. 63) and see if that gives a better understanding of how everything works.

http://robslink.com/SAS/book/GraphBeyond.pdf

I would recommend getting one thing working separately at a time (the ginside, the annotate, etc), and make sure you know how each piece works, before trying to get it all working together.

[Juno77]

Thanks. The proc ginside examples in your book were very helpful, the solution turned out to be simple. Now the problem is the map is not projecting correctly when all of the annotations are included.

[Juno77]

Which map dataset is the best to use in this case? I've been using maps.states, you used mapsgfk.us_states. What's the difference?

[GraphGuy]

The maps.* datasets are all older (and probably out of date a bit), and have the coordinates in radians.

The mapsgfk.* datasets are new/current, and have the coordinates in degrees.

Definitely use the newer mapsgfk datasets.

[Juno77]

Making progress using your example, but how do I add 2 more circles with increased distance from center zip? Thanks

[Juno77]

Re previous Q - Do I have to create 2 more circle datasets, 1 for each distance- that seems so messy.

[GraphGuy]

In this case, the circles are map areas, and a map can have multiple areas (each with a unique id) ... but, this becomes tricky when the circular areas are overlapping and/or concentric. Proc Gmap can 'draw' them fine, but when you're trying to do calculations (such as using proc ginside to determine which map area a point is inside of), then (I think) SAS procs assume the map areas are non-overlapping. Therefore, in this case, to be totally on the safe side, you'll probably need to at least handle the ginside calculations with 3 separate circle maps.

I haven't tried this particular scenario, and as they say "your mileage may vary"!