Si alguna vez os encontráis con el problema a resolver de saber si una coordenada, una dirección, una población, etc. se encuentra dentro de un polígono o de una área cualquiera que tienes definida por un listado de coordenadas esta pequeña macro hecha con SAS puede ayudarte.
¿Se encuentran estos puntos fuera o dentro del polígono?
Tenemos para este ejercicio los siguientes datos. Una tabla que contiene el polígono y otra que contiene el listado de coordenadas a evaluar.
data POLIGONO;
poligono="40.52,-4.32 40.55,-4.32 40.56,-4.29 40.53,-4.24 40.52,-4.26 40.52,-4.32";
run;
data COORDENADAS;
input coor_x coor_y;
datalines;
40.54548 -4.29331 /*dentro*/
40.56126 -4.34155 /*fuera*/
40.54 -4.32 /*se encuentra en el borde de una línea horizontal*/
40.52 -4.30 /*se encuentra en el borde de una línea vertical*/
40.49994 -4.28138 /*fuera*/
40.70045 -3.91955 /*fuera*/
40.52 -4.32 /*vértice*/
;
run;
Hay varias maneras de determinar si un punto se encuentra dentro o fuera de un polígono, utilizaremos una válida tanto para polígonos cóncavos como convexos, basada en el teorema de la curva de Jordan. De esta manera, si trazamos una línea cualquiera a partir de nuestro punto a evaluar, esta cortará el polígono un número impar de veces si el punto se encuentra dentro del mismo o 0 o un número par en caso contrario.
Lo primero que haremos será transponer el contenido del polígono que está representado en el formato siguiente (las coordenadas separados por coma y cada dupla separada por un espacio):
40.52,-4.32 40.55,-4.32 40.56,-4.29 40.53,-4.24 40.51,-4.26 40.52,-4.32
Para ello calculamos el número de vértices calculando el número de espacios con un compress y luego transponemos el contenido de forma que la salida es una tabla en la que cada registro contiene un segmento del polígono.
Así nuestro cálculo consistirá en saber si una línea horizontal trazada a partir de nuestro punto cruza una o varias veces. Utilizamos la siguiente macro:.
%macro poligono;
data POLIGONO;
set POLIGONO;
num_vertices = length(poligono) - length(compress(poligono));
do x=1 to num_vertices;
origen_x = input(scan(poligono,(x*2)-1,' ,'),8.2);
origen_y = input(scan(poligono,(x*2),' ,'),8.2);
destino_x = input(scan(poligono,((x+1)*2)-1,' ,'),8.2);
destino_y = input(scan(poligono,((x+1)*2),' ,'),8.2);
output;
end;
run;
%let dsid = %sysfunc(open(POLIGONO));
data SALIDA (drop=xinters);
retain id_alerta;
set COORDENADAS;
retain borde interseccion;
%let nobs = %sysfunc(attrn(&dsid,NOBS));
%let rc = %sysfunc(fetchobs(&dsid,1));
%do i=1 %to &nobs;
%let num_vertices = %sysfunc(getvarn(&dsid,%sysfunc(varnum(&dsid,num_vertices))));
%let x = %sysfunc(getvarn(&dsid,%sysfunc(varnum(&dsid,x))));
%let origen_x = %sysfunc(getvarn(&dsid,%sysfunc(varnum(&dsid,origen_x))));
%let origen_y = %sysfunc(getvarn(&dsid,%sysfunc(varnum(&dsid,origen_y))));
%let destino_x = %sysfunc(getvarn(&dsid,%sysfunc(varnum(&dsid,destino_x))));
%let destino_y = %sysfunc(getvarn(&dsid,%sysfunc(varnum(&dsid,destino_y))));
%let rc = %sysfunc(fetch(&dsid));
if &x=1 then do;
borde=0;
interseccion=0;
end;
if coor_x = &origen_x and coor_y = &origen_y then do;
resultado='VERTICE';
output;
end;
else if not(coor_x = &destino_x and coor_y = &destino_y) then do;
if &origen_y = &destino_y and
&origen_y = coor_y and
coor_x > min(&origen_x,&destino_x) and
coor_y <= min(&origen_y,&destino_y) and
coor_y <= max(&origen_y,&destino_y) and
coor_x <= max(&origen_x,&destino_x) and
&origen_y ne &destino_y then do;
xinters = (coor_y - &origen_y) * (&destino_x - &origen_x) / (&destino_y - &origen_y) + &origen_x;
if xinters = coor_x then borde=1;
if (&origen_x = &destino_x or coor_x <= xinters) then interseccion=interseccion+1;
end;
if &x = &num_vertices then do;
if mod(interseccion,2)=1 or borde=1 then do;
resultado='DENTRO';
output;
end;
else do;
resultado='FUERA';
borde=0;
interseccion=0;
output;
end;
end;
end;
%end;
run;
%let rc = %sysfunc(close(&dsid));
%mend;
%poligono;
Y el resultado que obtenemos es el siguiente:
Esta macro tiene la particularidad de que identifica correctamente cuando un punto se encuentra en un borde, pero en esos casos el indicador interseccion indicaría que se encuentran fuera del área del polígono. Hay que tener en cuenta esto y considerarlo en función de nuestros intereses. Yo, aquí, los estoy asumiendo como dentro.
Referencias:
Basado en el código de AssemblySys dataServices: https://assemblysys.com/es/algoritmo-punto-en-poligono
Programando con SAS 