%let name=umstead_park_development; filename odsout '.'; /* Related to: http://www.newsobserver.com/news/business/article178525016.html */ /* Here is the zoom-level, and the min/max x/y numbers of the desired OSM tiles ... (Alternatively, it would be a simple matter to programmatically determine the tiles needed to display a given list of lat/lon point data.) You can convert a lat/lon & zoom-level to the OSM x/y tile number using the following equation: --------- %let zoomlvl=13; data foo; lon_deg=-78.8013442; lat_deg=35.7983612; lat_rad=(atan(1)/45)*lat_deg; n=2**&zoomlvl; x = ((lon_deg + 180) / 360) * n; y = (1 - (log(tan(lat_rad) + (1/cos(lat_rad))) / constant('pi') )) / 2 * n; run; proc print data=foo; run; x, y = 2302.83 3222.55 */ %let zoomlvl=13; %let min_x=2301; %let max_x=2304; %let min_y=3219; /* top of map */ %let max_y=3222; /* bottom of map */ /* http://c.tile.openstreetmap.org/13/2301/3219.png http://tile.stamen.com/terrain-background/13/2301/3219.png https://maps.wikimedia.org/osm-intl/13/2301/3219.png %let tileserv=c.tile.openstreetmap.org; */ %let tileserv=a.basemaps.cartocdn.com/light_all; %include '../data_vault/proxy.sas'; /* I %include a file that sets the my_proxy macro variable. You won't be able to use my proxy, therefore set your own macro variable, which will be something like the following: %let my_proxy=http://yourproxy.com:80; */ /* Programmatically read in those map images, and save them to my hard drive so I can later annotate them (it would be nice if annotate could just refer to them without actually downloading them ... if anyone figures out how to do that, let me know! :) (This bit of code was adapted from some code Jim Metcalf was using to download a spreadsheet from the web. Note that your 'proxy' url will be different, for your own Intranet.)) */ %macro download_map(zoom, x, y); filename osm_pic url "http://&tileserv./&zoom./&x./&y..png" proxy="&my_proxy"; data _null_; n=-1; infile osm_pic recfm=s nbyte=n length=len _infile_=tmp; input; file "opnsta_&zoom._&x._&y..png" recfm=n; put tmp $varying32767. len; run; %mend; /* loop through all the desired map tile x/y numbers, and call the macro for each */ /* If you're re-running the same map several times, you probably want to comment this out, after the first time, since the png files are already downloaded. */ /* data _null_; do loop_x=&min_x to &max_x; do loop_y=&min_y to &max_y; call execute('%download_map('|| &zoomlvl ||','|| loop_x ||','|| loop_y ||');'); call execute('run;'); end; end; run; */ /* Create an annotate data set to 'stitch' the images together. */ data anno_tiles; length function style $8 html $200; xsys='2'; ysys='2'; when='b'; style='fit'; do loop_x=&min_x to &max_x; do loop_y=&min_y to &max_y; imgpath="opnsta_&zoomlvl._"||trim(left(loop_x))||"_"||trim(left(loop_y))||".png"; /* the y-coordinates might seem "reversed", but you will be "flipping" them later */ function='move'; x=loop_x; y=loop_y+1; output; function='image'; x=loop_x+1; y=loop_y; output; end; end; run; %let cout=gray33; /* outline color of land area */ /* %let folder=\\l7a695\public\wake\2017\; */ %let folder=d:\public\wake\2017\; libname robsdata "&folder"; data mymap; set robsdata.umstead_area; run; /* Use gproject to 'clip' out a certain area of the map, so I don't have to work with such a large dataset. (this is optional) */ /* proc gproject data=mymap out=mymap latlong eastlong degrees project=none latmax=35.905 latmin=35.810 longmin=-78.83 longmax=-78.72 ; id pin_num; run; */ /* Instead of projecting the map with Proc Gproject ... Convert the long & lat degrees to x & y in the OSM tile-based coordinate system, using equations from: http://wiki.openstreetmap.org/wiki/Tilenames SAS doesn't have a sec() function, so I used 1/cos(value) */ data mymap; set mymap; lon_rad=(atan(1)/45)*long; lat_rad=(atan(1)/45)*lat; n=2**&zoomlvl; x = ((long + 180) / 360) * n; y = (1 - (log(tan(lat_rad) + (1/cos(lat_rad))) / constant('pi') )) / 2 * n; run; data anno_circles_labels; length label_position $1 label_text $100; infile datalines dlm=','; input lat long pie_size label_angle label_position label_text; datalines; 35.8408933,-78.7757316,3.8,70,3,Proposed New Quarry 35.8449202,-78.7947162,2.5,250,d,Proposed Luxury Hotel / Offices 35.8483811,-78.7839606,4.8,90,b,Proposed Parking ; run; data anno_circles_labels; set anno_circles_labels; lon_rad=(atan(1)/45)*long; lat_rad=(atan(1)/45)*lat; n=2**&zoomlvl; x = ((long + 180) / 360) * n; y = (1 - (log(tan(lat_rad) + (1/cos(lat_rad))) / constant('pi') )) / 2 * n; run; /* Since openstreetmaps starts at top/left, and sas starts at bottom/left, "flip" the y-coordinate. */ data mymap; set mymap; y=y*-1; run; data anno_tiles; set anno_tiles; y=y*-1; run; data anno_circles_labels; set anno_circles_labels; y=y*-1; run; /* These are the circled areas, and their labels */ data anno_circles_labels; set anno_circles_labels; length function $8 style $35 color $8; xsys='2'; ysys='2'; hsys='3'; when='a'; /* the circle */ function='pie'; width=2; style='pempty'; rotate=360; color='blue'; size=pie_size; output; /* the label outside the circle */ function='pie'; style='pempty'; line=0; size=size*1.1; angle=label_angle; rotate=0; output; function='piexy'; size=1; output; function='cntl2txt'; output; function='label'; style=''; angle=0; rotate=0; size=.; x=.; y=.; color='blue'; size=1.6; position=label_position; text=trim(left(label_text)); output; run; data anno_text; length function $8 color $20 style $35 text $100; xsys='3'; ysys='3'; hsys='3'; when='a'; function='label'; color=''; style='albany amt/bold'; position='5'; cbox='cxf5f5f3'; x=50; y=96; size=5; text="Umstead Park Development"; output; run; /* Combine all the annotate datasets */ data anno_all; set anno_text anno_circles_labels; run; data mydata; set robsdata.umstead_area_data; length my_html $500 colorval $50; if index(propdesc,"WILLIAM UMSTEAD STATE PARK")^=0 then colorval="Umstead Park"; /* PINs provided by Jean Spooner, umsteadcoalition.org */ if pin_num in ('0776275726' '0775364221' '0775542715' '0767993962' '0768914024' '0778007683' '0777452678' '0777344265' '0777424290' '0777700300' '0777704352' '0785187931' '0785187931') then colorval="Umstead Park"; if index(propdesc,"RALEIGH-DURHAM INTN'L AIRPORT")^=0 then colorval="Airport"; if index(propdesc,"PROP OF RALEIGH DURHAM AIRPORT")^=0 then colorval="Airport"; if index(owner,"AIRPORT AUTHORITY")^=0 then colorval="Airport"; if index(owner,"WAKE STONE CORPORATION")^=0 and index(propdesc,"YOUNG LAND")^=0 then colorval="Quarry"; if index(owner,"SAS INSTITUTE INC")^=0 then colorval="SAS"; if index(propdesc,"CRABTREE CREEK WATERSHED PROJECT SITE 23")^=0 then colorval="Crabtree Park"; if colorval^='' then do; my_html='title='||quote( 'Parcel: '||trim(left(pin_num))||'0d'x|| 'Owner: '||trim(left(owner))||'0d'x|| 'Propdesc: '||trim(left(propdesc)))|| ' href='||quote('http://services.wakegov.com/realestate/Account.asp?id='||trim(left(reid))); output; end; run; goptions device=png; goptions xpixels=850 ypixels=900; goptions noborder; ODS LISTING CLOSE; ODS HTML path=odsout body="&name..htm" (title="Umstead Park, NC development") style=htmlblue; goptions gunit=pct htext=1.8 ftext="albany amt/bold"; goptions ctext=gray33; legend1 label=none mode=protect position=(top right inside) across=1 shape=bar(.15in,.15in) offset=(-4,0) cframe=cxf5f5f3 cborder=gray77 order=('Airport' 'Umstead Park' 'Quarry' 'Crabtree Park' 'SAS'); /* Alpha-transparent colors for the land parcels */ pattern1 v=s c=Aff7f0033; pattern2 v=s c=Affff3333; pattern3 v=s c=Ae41a1c77; pattern4 v=s c=A984ea377; pattern5 v=s c=A00ff0033; title1 h=12 " "; /* I only want certain parcels to show up, so I don't use gmap's "all" option */ proc gmap data=mydata map=mymap anno=anno_tiles; note font='albany amt' link='http://www.cartocdn.com/copyright' move=(20,0.5) height=8pt 'a9a0'x color=blue underlin=1 "OpenStreetMap" c=gray33 underlin=0 " contributors"; id pin_num; choro colorval / legend=legend1 coutline=gray55 cempty=graydd anno=anno_all html=my_html des='' name="&name"; run; quit; ODS HTML CLOSE; ODS LISTING;