%let gpath='.';
%let dpi=200;
ods html close;
ods listing style=htmlblue gpath=&gpath image_dpi=&dpi;
/*--Define some random polygons--*/
data polygons;
PI=constant("PI");
xmin=1e6; ymin=1e6;
npoly=3;
do pid=1 to npoly;
r=10+10*ranuni(3);
nodes=int(3+3*ranuni(2));
xoff=ifn(xmax ne ., xmax, 0)+r;
yoff=ifn(ymax ne ., ymax, 0)+r;
step=360/(nodes);
nid=0;
do deg=0 to 359 by step;
nid=nid+1;
x=xoff+(r*(1+ranuni(2))*cos(deg*pi/180));
y=yoff+(r*(1+ranuni(2))*sin(deg*pi/180));
xmax=max(x, xmax);
ymax=max(y, ymax);
xmin=min(x, xmin);
ymin=min(y, ymin);
output;
end;
end;
/*--Hand build a "L" polygon--*/
pid=npoly+1;
nodes=6;
id=1;
nid=1; x=50; y=30; xmin=min(x, xmin); ymin=min(y, ymin); output;
nid=2; x=50; y= 0; xmin=min(x, xmin); ymin=min(y, ymin); output;
nid=3; x=70; y= 0; xmin=min(x, xmin); ymin=min(y, ymin); output;
nid=4; x=70; y=10; xmin=min(x, xmin); ymin=min(y, ymin); output;
nid=5; x=60; y=10; xmin=min(x, xmin); ymin=min(y, ymin); output;
nid=6; x=60; y=30; xmin=min(x, xmin); ymin=min(y, ymin); output;
/*--Hand build a Right Triangle--*/
pid+1;
nodes=3;
id=2;
nid=1; x=80; y=50; xmin=min(x, xmin); ymin=min(y, ymin); output;
nid=2; x=80; y=10; xmin=min(x, xmin); ymin=min(y, ymin); output;
nid=3; x=120; y=10; xmin=min(x, xmin); ymin=min(y, ymin); output;
call symput ("XMin", xmin);
call symput ("YMin", ymin);
run;
/*%put &xmin &ymin;*/
/*ods html;*/
/*proc print data=polygons;*/
/*var pid nodes nid x y;*/
/*run;*/
/*ods html close;*/
/*--Plot polygons--*/
ods graphics / reset width=4in height=3in imagename='polyAreaSG';
title 'Polygons';
proc sgplot data=polygons noborder noautolegend;
polygon id=pid x=x y=y / fill group=pid dataskin=sheen;
xaxis display=none min=&xmin;
yaxis display=none;
run;
/*--Define Macro to compute polygonal area and center of area--*/
%macro polyarea (ds=, xmin=, ymin=, out=, Id=, X=, Y=);
data &out;
retain xp yp ax ay xm ym xfirst yfirst n;
set &ds(keep=&id &x &y);
by &id;
area=.; xc=.; yc=.; xo=.; yo=.;
/*--Process first node--*/
if first.&id then do;
n=1; /*--Node id --*/
ax=0; /*--Sum area for x--*/
ay=0; /*--Sum area for y--*/
dax=0; /*--Area of segment for x--*/
day=0; /*--Area of segment for y--*/
xm=0; /*--Moment of segment area about minx--*/
ym=0; /*--Moment of segment area about miny--*/
xfirst=&x; yfirst=&y;
output;
end;
/*--Process each polygon segment--*/
else do;
/*--Area of segment based on X parallelogram, Moment about XMin & sum of X areas--*/
dax=(&x-&xmin+xp-&xmin)*(&y-yp)/2;
a=&x-&xmin;
b=xp-&xmin;
v=a+b;
h=&y-yp;
dx=0; y1=0;
/*--Compute height of centroid from base of trapezoid--*/
if h ne 0 and v ne 0 then do;
y1=h*(2*a+b)/(3*(a+b));
dx=a/2+(b-a)*(h-y1)/(2*h);
end;
/*--Compute moment for atea about xmin--*/
xm=xm+dax*dx;
ax=ax+dax;
/*--Compute moment for atea about ymin--*/
dy=yp+y1-&ymin;
ym=ym+dax*dy;
n+1;
output;
end;
xp=&x; yp=&y;
/*--Process final polygon segment--*/
if last.&id then do;
dax=(xfirst-&xmin+&x-&xmin)*(yfirst-&y)/2;
a=xfirst-&xmin;
b=&x-&xmin;
v=a+b;
h=yfirst-&y;
dx=0; y1=0;
/*--Compute height of centroid from base of trapezoid--*/
if h ne 0 and v ne 0 then do;
y1=h*(2*a+b)/(3*(a+b));
dx=a/2+(b-a)*(h-y1)/(2*h);
end;
/*--Compute moment for atea about xmin--*/
xm=xm+dax*dx;
ax=ax+dax;
/*--Compute moment for atea about ymin--*/
dy=&y+y1-&ymin;
ym=ym+dax*dy;
area=abs(ax);
/*--Compute (x, y) of centroid--*/
xc=abs(xm)/ area + &xmin;
yc=abs(ym)/ area + &ymin;
xp=&x; yp=&y;
&x=xfirst; &y=yfirst;
output;
end;
run;
%mend;
/*--Invoce macro--*/
%polyarea (ds=polygons, xmin=&xmin, ymin=&ymin, out=area, id=pid, x=x, y=y);
/*proc print data=area;*/
/*var pid xp yp x y ax ay area dx dy dax day xm ym xc yc y1;*/
/*run;*/
/*--Add x and y values for labels--*/
data arealbl;
set area;
xl=xc; yl=yc+4;
run;
/*--Draw polygons with labeled area and CG--*/
ods graphics / reset width=4in height=3in imagename='polyAreaCentroidSG';
title 'Polygons with Area and Centroid';
proc sgplot data=arealbl noborder noautolegend;
format area f5.0;
polygon id=pid x=x y=y / fill group=pid dataskin=sheen;
text x=xl y=yl text=area / group=pid strip;
scatter x=xc y=yc / markerattrs=(symbol=circle size=3);
xaxis display=none;
yaxis display=none;
run;
title;
proc template;
define statgraph polygons;
begingraph;
entrytitle 'Polygons with Area and Centroid';
layout overlayequated / walldisplay=none cycleattrs=false
xaxisopts=(display=none) yaxisopts=(display=none);
polygonplot id=pid x=x y=y / display=(fill) group=pid dataskin=sheen;
textplot x=xl y=yl text=eval("A=" || put(area,4.0)) / group=id strip=true;
scatterplot x=xc y=yc / markerattrs=(symbol=circle size=3) group=pid;
endlayout;
endgraph;
end;
run;
ods graphics / reset width=4in height=3in imagename='polyAreaCentroidEqGTL';
proc sgrender data=arealbl template=polygons;
format area 4.0;
run;