BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.
genemroz
Quartz | Level 8

Recently ballardw helped me with a question about how to test whether points lie within a polygon using PROC GINSIDE.

In the code below, I have modified the example code he provided to include a third polygon that overlaps the two he created for his example. Now, two points, #4 and #6, each lie within two polygons. Point #4 lies within polygons 1 and 3. Point # 5 lies within polygons 2 and 3. PROC GINISDE reports that these two points lies within polygon 1 and 2, respectively, but does not report that they are also within polygon 3.

How can I modify the code so that PROC GINISDE correctly reports all polygons occupied by a point?  If possible, I would like to retain all the polygons within a single map file.

Thanks in advance for any help,

Gene

data examplemap;
   input id x y;
datalines;
1   0 0
1   0 10
1   5 15
1   8 10
1   8 0
1   0 0
2   12 12
2   12 20
2   18 20
2   18 12
2   12 12
3   3  9
3   3  19
3   16  19
3   16 9
3   3  9
;


data exampledata;
   input otherid x y;
datalines;
1  0  .5
2  1  4
3  0  20
4  4  12
5  10 10
6  13 14
7  20 1
;

proc ginside map =examplemap
             data=exampledata
             out =result
             includeborder
;
   id id;
run;
1 ACCEPTED SOLUTION

Accepted Solutions
GraphGuy
Meteorite | Level 14

This is a bit of a 'grey' area...

 

Although Proc Gmap lets you draw overlapping polygons, I think Proc Ginside relies on the underlying assumption that all the map polygons are "topologically correct" (no overlapping polygons, and all common borders have common points in both polygons, etc).

 

To get the results you want, I think you'd have to run a separate Proc Ginside for each polygon in the map (if you have overlapping polygons).

 

View solution in original post

4 REPLIES 4
GraphGuy
Meteorite | Level 14

This is a bit of a 'grey' area...

 

Although Proc Gmap lets you draw overlapping polygons, I think Proc Ginside relies on the underlying assumption that all the map polygons are "topologically correct" (no overlapping polygons, and all common borders have common points in both polygons, etc).

 

To get the results you want, I think you'd have to run a separate Proc Ginside for each polygon in the map (if you have overlapping polygons).

 

ballardw
Super User

This creates some additional variables in your exampledata to use it as a map annotate data set and then shows the values on a map.

data exampledata;
   input otherid x y ;
   text=put(otherid,1.); 
   color='black'; 
   size=2; 
   position='5'; 
retain xsys ysys '2' hsys '3' when 'a'   ;
datalines;
1  0  .5
2  1  4
3  0  20
4  4  12
5  10 10
6  13 14
7  20 1
;

proc gmap map=examplemap data=examplemap annotate=exampledata;
 id  id ;
 choro id  ;
run;
quit;
  

Gmap doesn't do transparency but when I look at the plotted values of your example data I do not see anything that looks like the three areas even overlap much less that any of your points are in all three areas.

 

I suggest that you get a piece of graph paper and trace out the coordinates of your map areas (carefully) so you can see them clearly and then plot the points in question.

GraphGuy
Meteorite | Level 14

One correction/addition to what ballardw said - gmap does support transparency. Just specify an alpha-transparent color on the pattern statement (see pattern3 in this example, for the green polygon):

 

data examplemap;
input id x y;
datalines;
1 0 0
1 0 10
1 5 15
1 8 10
1 8 0
1 0 0
2 12 12
2 12 20
2 18 20
2 18 12
2 12 12
3 3 9
3 3 19
3 16 19
3 16 9
3 3 9
;

 

data exampledata;
input otherid x y;
datalines;
1 0 .5
2 1 4
3 0 20
4 4 12
5 10 10
6 13 14
7 20 1
;

 

data anno_dots; set exampledata;
xsys='2'; ysys='2'; hsys='3'; when='a';
function='pie'; style='psolid'; rotate=360; size=1.0;
run;

 

pattern1 v=s c=red;
pattern2 v=s c=blue;
pattern3 v=s c=A00ff0055;

 

proc gmap data=examplemap map=examplemap anno=anno_dots;
id id;
choro id;
run;

 

proc ginside map=examplemap data=exampledata out=result includeborder;
id id;
run;

 

proc print data=result; run;

 

inside.png

 

inside_table.png

genemroz
Quartz | Level 8

Robert, 

Thanks for this response.  I was certain I had set up the example so that points 4 and 6 were in both polygons as the output of your code shows.

 

If I pursue this topic of identifiying points within overlapping polygons, I'll have to figure out how to loop through multiple calls to PROC GINSIDE to access many (hundreds) of polygons. It may be more trouble than its worth.

 

Thanks again,

 

Gene

 

sas-innovate-2024.png

 

Time is running out to save with the early bird rate. Register by Friday, March 1 for just $695 - $100 off the standard rate.

 

Check out the agenda and get ready for a jam-packed event featuring workshops, super demos, breakout sessions, roundtables, inspiring keynotes and incredible networking events. 

 

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.

Click image to register for webinarClick image to register for webinar

Classroom Training Available!

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

View all other training opportunities.

Discussion stats
  • 4 replies
  • 667 views
  • 1 like
  • 3 in conversation