/* SAS program to accompany the article "Offset regions: Find all points within a specified distance from a polygon" by Rick Wicklin, published 16JUL2018 on The DO Loop blog: https://blogs.sas.com/content/iml/2018/07/16/points-distance-from-polygon.html This program shows how to construct the "offset region" for a simple polygon. The offset region of radius r is the union of - The polygon itself - Circles of radius r centered at each vertex - Rectangles of width r positions outside the polygon along each edge */ /* Before starting the SAS/IML program, generate a polygon for the state of Texas in standardized coordinates. */ data Texas1; set maps.US; where StateCode = "TX" and Segment=1; by State Segment; run; * Create density variable that can be used to reduce the data set; proc greduce data=Texas1 out=Texas2; id segment; run; data Texas; set Texas2(keep=segment x y density); where (density<=4); Obs = _N_; run; /************************************************************************/ /* Define functions to compute and plot the offset region for a polygon */ /************************************************************************/ title; ods graphics / reset width=350px height=350px subpixel=on; proc iml; /* Construct a union of polygons. Place a circle at every vertex and a rectangle along each edge. */ start OffsetUnion(Poly, r); M = {0 -1, 1 0}; /* rotation matrix by 90 degrees */ t = 2*constant('pi')*T( do(0, 1, 0.01) ); /* parameterization of [0, 2*pi] */ k = nrow(t); P = Poly; N = nrow(P); P = Poly // Poly[1,]; /* for convenience, append first point to end */ do i = 1 to N; id = j(k, 1, 2*i-1); /* ID for i_th circle (constant vector) */ circle = P[i,] + (r*cos(t) || r*sin(t)); /* circle of radius r centered at i_th vertex */ G = G // (id || circle); /* append polygon approx of circle */ v = P[i+1, ] - P[i, ]; /* vector from i_th to (i+1)_th vertex */ nv = r / norm(v) * v*M ; /* scaled normal to v */ p1 = P[i,] + nv; /* offset of i_th vertex */ p2 = P[i+1,] + nv; /* offset of (i+1)_th vertex */ id = j(4, 1, 2*i); /* ID for i_th rectangle */ G = G // (id || (p1 // p2 // P[i+1,] // P[i,])); /* append rectangle */ end; return G; finish; start PlotOffset(G, P, t2); create Offset from G[c={ID x y}]; append from G; close; PP = j(nrow(P), 1, 0) || P; create Poly from PP[c={ID2 x y}]; append from PP; close; submit t2; data All; set Offset Poly; run; title "Union of Circles and Rectangles"; title2 "&t2"; proc sgplot data=All noautolegend; polygon id=ID2 x=x y=y / fill fillattrs=GraphData2; polygon id=ID x=x y=y; run; title "Offset Region"; title2 "&t2"; proc sgplot data=All noautolegend; polygon id=ID x=x y=y / fill; polygon id=ID2 x=x y=y / fill outline fillattrs=GraphData2; run; endsubmit; finish; P = {1 0, /* convex polygon */ 1 1, 0 2, -1 1, -1 0}; r = 0.5; G = OffsetUnion(P, r); run PlotOffset(G, P, "r = 0.5"); P = {2 0, /* nonconvex C-shaped with two prongs that are near */ 2 2, 0.8 2.5, 1.5 1.5, 1 1, 0 1, -0.5 1.5, 0.4 2.5, -1 2, -1 0}; r = 0.4; G = OffsetUnion(P, r); run PlotOffset(G, P, "r = 0.4"); use Texas; read all var {x y} into P; close; r = max( distance(P) ); print r; G = OffsetUnion(P, r); run PlotOffset(G, P, "r = 0.2036");