/*-----------------------------------------------------------------* | Name: gfrance.sas | | Title: Create maps of France ~ 1830 corresponding to Guerry | | from maps.france | | ----------------------------------------------------------------| | Procs: sort print format gremove greduce contents copy gmap | | Macdefs: fixdepts copyem gexplode shift | | Macros: mixcase centroid | | Datasets: gfrance gfrance2 gregion goutline mymaps.xfrance | | mymaps.xfrance2 mymaps.xoutline names | | ----------------------------------------------------------------| | Author: Michael Friendly | | Created: 19 Jun 2006 10:59:00 | *-----------------------------------------------------------------*/ /* title: Create maps of France ~ 1830 corresponding to Guerry from maps.france */ libname mymaps '~/sasuser/maps'; *-- Fix the map of France to conform to Guerry: - adjust the 97 current departments to correspond to the 86 in 1830 - delete those not part of France then ; %macro fixdepts; /* Corse was one dept - merge these to one area, new ID */ if id in (201, 202) then Dept=200; /* Seine et Oise (78) was cut into Essonne (91), Val d'Oise (95) and Yvelines (78) */ else if id in (91, 95) then Dept=78; /* Seine (75) now split into Hauts-de-Seine (92), Seine-Saint-Denis (93) et Val-de-Marne (94)*/ else if id in (92, 93, 94) then Dept=75; /* departments not part of France in 1830 */ else if id in ( 6, /* Alpes-Maritimes */ 73,74, /* Savoie, Haute-Savoie */ 90) /* Territore-de-Belfort */ then delete; else Dept=id; %mend; data gfrtemp; set maps.france; %fixdepts; run; *-- remove internal boundaries based on merged DEPT; proc sort data=gfrtemp; by dept; proc gremove data=gfrtemp out=gfrance; by dept; id id; run; *-- Add density value to data set (0-6); proc greduce data=gfrance out=gfrance; id dept; *-- shift Corsica so it doesn't extend the BBox too much; %let shift=25; data gfrance(label='Map of France, 1830'); set gfrance; if dept=200 then do; x = x-&shift; y = y+&shift; end; label Dept = 'Department Id number' /* Department = 'Department name' */ Density = 'Map density value' segment = 'Map polygon segment number' ; *-- Fix the map labels for France to conform to Guerry: ; data gfrance2(label='France department labels and centroid'); set maps.france2; drop _MAP_GEOMETRY_ COUNTRY REGION; %mixcase(idname,dlim="-'"); *mixcase(region,dlim="-'"); %fixdepts; proc sort data=gfrance2; by dept; /* Read Guerry data - collapse the 22 current regions to the 5 considered by Guerry (N S E W Central) */ *include 'guerry_a2.sas'; %include data(guerry); data gregion; * set guerry_a2; set guerry; keep dept Region Department; proc sort data=gregion; by dept; *-- Find centroids of departments; %centroid(gfrance, dept, out=gcenter, xc=x, yc=y); data gfrance2(label='France department labels and centroid'); merge gfrance2 gregion gcenter; by dept; if first.dept; if dept=200 then do; idname='Corse'; id=200; end; * drop regionid; *proc contents data=gfrance2; proc print data=gfrance2; id dept; *-- use a data step to convert the data set gfrance into an Annotate data set goutline - specify thick lines ; data fmt / view=fmt; retain fmtname 'dept_'; set guerry (rename=(dept=start region=label)); keep fmtname start label; run; proc format cntlin=fmt; run; data goutline; set gfrance; where dept ne 200; region = put(dept,dept_.); run; proc sort data=goutline; by region; run; proc gremove data=goutline out=goutline; by region; id dept; run; *-- Add density value to data set (0-6); proc greduce data=goutline out=goutline; id region; label Density = 'Map density value' ; data goutline(label='Region outline annotate dataset'); length function color $8; set goutline; by region segment; retain xsys ysys '2' hsys '3' size 2 color 'black' when 'a'; if first.segment or (lag(x)=. and lag(y)=.) then function='POLY'; else function='POLYCONT'; if x and y then output; drop segment; run; proc contents data=goutline; %macro copyem; proc copy in=work out=mymaps memtype=data; select gfrance gfrance2 goutline; run; %mend; %copyem; %macro gexplode(offset=0); *-- Explode the regions of France for visibility; %macro shift; select(region); when('N') do; y = y+&offset; end; when('S') do; y = y-&offset; end; when('E') do; x = x+&offset; end; when('W') do; x = x-&offset; end; otherwise; end; %mend; %if &offset>0 %then %do; data mymaps.xfrance(label="Map of France, 1830, regions exploded by &offset"); merge mymaps.gfrance mymaps.gfrance2(keep=dept region); by dept; %shift; data mymaps.xfrance2(label="ID names and centroids, regions exploded by &offset"); set mymaps.gfrance2; %shift; run; data mymaps.xoutline(label="Region outline annotate dataset, regions exploded by &offset"); set mymaps.goutline; %shift; run; %end; %mend; %gexplode(offset=8); %let prefix=x; %let lib=mymaps; data names; set &lib..&prefix.france2; length function $8 text $20; xsys ='2'; ysys='2'; when='A'; function='label'; text=idname; size=0.7; position='E'; output; text = put(dept, 3.); size=1; position='B'; output; data names; set names &lib..&prefix.outline; /* Test mymaps.xfrance */ %include goptions; pattern1 v=solid c=cx99CC66; * C; pattern2 v=solid c=cx666699; * E; pattern3 v=solid c=cxCC9966; * N; pattern4 v=solid c=cxCC6699; * S; pattern5 v=solid c=cx9966CC; * W; proc gmap map=&lib..&prefix.france (where=(density<=4)) data=&lib..&prefix.france2; id dept; /* choro idname / discrete coutline=graycc name='gfrance' nolegend */ choro region / discrete coutline=graycc name="gfrance" nolegend anno=names ;