BookmarkSubscribeRSS Feed
☑ This topic is solved. Need further help from the community? Please sign in and ask a new question.
Juno77
Calcite | Level 5

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;

 

1 ACCEPTED SOLUTION

Accepted Solutions
GraphGuy
Meteorite | Level 14

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.

View solution in original post

9 REPLIES 9
GraphGuy
Meteorite | Level 14

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

 

customers_within_circle.png

Juno77
Calcite | Level 5

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
Meteorite | Level 14

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
Calcite | Level 5
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
Calcite | Level 5
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
Meteorite | Level 14

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
Calcite | Level 5
Making progress using your example, but how do I add 2 more circles with increased distance from center zip? Thanks
Juno77
Calcite | Level 5
Re previous Q - Do I have to create 2 more circle datasets, 1 for each distance- that seems so messy.
GraphGuy
Meteorite | Level 14

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"!

 

SAS INNOVATE 2024

Innovate_SAS_Blue.png

Registration is open! SAS is returning to Vegas for an AI and analytics experience like no other! Whether you're an executive, manager, end user or SAS partner, SAS Innovate is designed for everyone on your team. Register for just $495 by 12/31/2023.

If you are interested in speaking, there is still time to submit a session idea. More details are posted on the website. 

Register now!

How to Concatenate Values

Learn how use the CAT functions in SAS to join values from multiple variables into a single value.

Find more tutorials on the SAS Users YouTube channel.

Get the $99 certification deal.jpg

 

 

Back in the Classroom!

Select SAS Training centers are offering in-person courses. View upcoming courses for:

View all other training opportunities.

Discussion stats
  • 9 replies
  • 400 views
  • 3 likes
  • 2 in conversation