Last week I posted a question entitled "Heatmap with geographic features?" Robert Allison came up with a suggested path to a solution that I've attached below. The code for this solution is also below.
Now I'd like to build on that solution to create two additional maps:
1.) A similar heat map but limited to only a few counties within the state boundaries.
2.) A similar heat map but one that extends beyond the New Mexico state boundaries to include county lines (and a few cities) in adjoining states (AZ, CO, OK, TX) but only extending into these states by about 1 degree of latitude and longitude.
Basically both questions are about how to define the area covered by county polygons according to user-selected latitude and longitude limits. I tried using GPROJECT to "clip" the polygons but this didn't seem to work--maybe I wasn't using it correctly.
Thanks in advance for any help.
Gene
data my_map; set mapsgfk.us_counties (where=(density<=3 and statecode='NM'));
length id_plus_segment $50;
id_plus_segment=trim(left(id))||'_'||trim(left(segment));
run;
data city_data; set mapsgfk.uscity (where=(statecode='NM' and pop_type^='under10k'
or (statecode='NM' and city='Socorro') or (statecode='NM' and city='Quemado') or
(statecode='NM' and city='Magdalena') or (statecode='NM' and city='Grants')
or (statecode='NM' and city="Bernalillo")));
if city= "North Valley" then delete;
if city="South Valley" then delete;
city_long=long;
city_lat=lat;
citylabel=city;
run;
data all_csrnmdata; set my_map city_data dec20.csrnm;
run;
title h=12pt "Meteors Observed by NMMA: December 2020";
proc sgplot data=all_csrnmdata noborder noautolegend;
polygon x=long y=lat id=id_plus_segment / tip=none
nofill outline lineattrs=(color=black);
scatter x=city_long y=city_lat / markerattrs=(color=blue)
datalabel=citylabel datalabelpos=bottom datalabelattrs=(color=blue);
heatmap x=lng1 y=lat1 / name="heatmap" nxbins=50 nybins=50
outline outlineattrs=(color=graycc thickness=0) transparency=.6
colormodel=(cxfeebe2 cxfbb4b9 cxf768a1 cxc51b8a cx7a0177);
gradlegend /title='Count';
xaxis display= (nolabel);
yaxis display= (nolabel);
run;
You will need to provide some example data for dec20.csrnm if you want us to test code.
See if the example of using MIN= and MAX= on the XAXIS and YAXIS statements is doing similar to what you want.
title h=12pt "Plot without heatmap"; proc sgplot data=all_csrnmdata noborder noautolegend; polygon x=long y=lat id=id_plus_segment / tip=none nofill outline lineattrs=(color=black); scatter x=city_long y=city_lat / markerattrs=(color=blue) datalabel=citylabel datalabelpos=bottom datalabelattrs=(color=blue); /*heatmap x=lng1 y=lat1 / name="heatmap" nxbins=50 nybins=50*/ /*outline outlineattrs=(color=graycc thickness=0) transparency=.6*/ /* colormodel=(cxfeebe2 cxfbb4b9 cxf768a1 cxc51b8a cx7a0177);*/ gradlegend /title='Count'; xaxis display= (nolabel) min=-106.2 max=-104.7; yaxis display= (nolabel) min=33 max=36; run;title;
You would select specific counties with a where like
Where county in ( 3 , 4 , 5)
which I realize probably won't make geographic sense.
Or make sure that the boundary data has the County_name set for every record and then list the desired county names
Where county_name in ("Bernalillo" "Santa Fe" "Soccoro")
These could be a statement or a data set option.
You will need to provide some example data for dec20.csrnm if you want us to test code.
See if the example of using MIN= and MAX= on the XAXIS and YAXIS statements is doing similar to what you want.
title h=12pt "Plot without heatmap"; proc sgplot data=all_csrnmdata noborder noautolegend; polygon x=long y=lat id=id_plus_segment / tip=none nofill outline lineattrs=(color=black); scatter x=city_long y=city_lat / markerattrs=(color=blue) datalabel=citylabel datalabelpos=bottom datalabelattrs=(color=blue); /*heatmap x=lng1 y=lat1 / name="heatmap" nxbins=50 nybins=50*/ /*outline outlineattrs=(color=graycc thickness=0) transparency=.6*/ /* colormodel=(cxfeebe2 cxfbb4b9 cxf768a1 cxc51b8a cx7a0177);*/ gradlegend /title='Count'; xaxis display= (nolabel) min=-106.2 max=-104.7; yaxis display= (nolabel) min=33 max=36; run;title;
You would select specific counties with a where like
Where county in ( 3 , 4 , 5)
which I realize probably won't make geographic sense.
Or make sure that the boundary data has the County_name set for every record and then list the desired county names
Where county_name in ("Bernalillo" "Santa Fe" "Soccoro")
These could be a statement or a data set option.
I think you could modify your map data using GPROJECT with option PROJECT=NONE and options LATMIN=, LATMAX=, LONGMIN=, and LONGMAX=
Esteemed Advisers:
I'm re-opening this thread because I've encountered a new difficulty.
Robert Allison and Ballardw recently helped me create a heatmap that overlays a state map using SGPLOT. And it works fine until I tried to add some labels to the markers that denote location of selected cities. Then the aspect ratio gets wonky (a technical term for "compressed in the horizontal direction"). I thought the solution was to use the STRIP function to remove any leading or trailing spaces from the city variable but that didn't seem to help.
Attached to this post is the output of the code below. It shows the correct and desired map without datalabels for city names and below that is what results when datalabels for city names are included. I've omitted the code for generating the heatmap because it's not relevant to the issue of how/why adding datalabels impacts the aspect ratio.
Thanks in advance for any guidance you can provide.
Gene
data my_map; set mapsgfk.us_counties (where=(density<=3 and statecode='NM'));
length id_plus_segment $50;
id_plus_segment=trim(left(id))||'_'||trim(left(segment));
run;
data city_data; set mapsgfk.uscity (where=(statecode='NM' and pop_type^='under10k'
or (statecode='NM' and city='Socorro') or (statecode='NM' and city='Quemado') or
(statecode='NM' and city='Magdalena') or (statecode='NM' and city='Grants') or
(statecode='NM' and city='Truth or Consequences')or (statecode='NM' and city="Bernalillo")or
(statecode='NM' and city="Taos") or (statecode='NM' and city="Tucumcari")));
if city= "North Valley" then delete;
if city="South Valley" then delete;
if city='Sunland Park' then delete;
city_long=long;
city_lat=lat;
citylabel=strip(city);
run;
data fovbox;
input Corners $1-3 fovx fovy;
datalines;
fov -108.820 33.30
fov -108.820 36.90
fov -104.430 36.90
fov -104.430 33.30
fov -108.820 33.30
;
run;
data all_csrnmdata; set my_map city_data fovbox;
run;
title h=12pt "City Labels Excluded";
proc sgplot data=all_csrnmdata noautolegend aspect=1.221;
polygon x=long y=lat id=id_plus_segment / tip=none
nofill outline lineattrs=(color=black);
polygon x=fovx y=fovy id=corners/outline lineattrs=(color=red thickness=1);
scatter x=city_long y=city_lat / markerattrs=(color=blue);
refline -106.6293 / axis=x lineattrs=(thickness=1 color=blue);
refline 35.1055 / axis=y lineattrs=(thickness=1 color=blue);
gradlegend /title='Count';
xaxis display= (nolabel);
yaxis display= (nolabel);
run;
title h=12pt "City Labels Included";
proc sgplot data=all_csrnmdata noautolegend aspect=1.221;
polygon x=long y=lat id=id_plus_segment / tip=none
nofill outline lineattrs=(color=black);
polygon x=fovx y=fovy id=corners/outline lineattrs=(color=red thickness=1);
scatter x=city_long y=city_lat / markerattrs=(color=blue)
datalabel=citylabel datalabelpos=top datalabelattrs=(color=blue);
refline -106.6293 / axis=x lineattrs=(thickness=1 color=blue);
refline 35.1055 / axis=y lineattrs=(thickness=1 color=blue);
gradlegend /title='Count';
xaxis display= (nolabel);
yaxis display= (nolabel);
run;
Registration is now open for SAS Innovate 2025 , our biggest and most exciting global event of the year! Join us in Orlando, FL, May 6-9.
Sign up by Dec. 31 to get the 2024 rate of just $495.
Register now!
Learn how use the CAT functions in SAS to join values from multiple variables into a single value.
Find more tutorials on the SAS Users YouTube channel.
Ready to level-up your skills? Choose your own adventure.