Can PROC GINSIDE test for points within overlapping polygons?

[genemroz]

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;

[GraphGuy · Accepted solution]

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]

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]

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;

[genemroz]

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