/* Download North Carolina shapefiles from here:
   http://data.nconemap.com/geoportal/catalog/main/home.page
   Search for Major Hydrography.
   Run the program below to import the shapefile and use SGMap   */

ods _all_ close;
ods html file="capefear.html";

/* import the map data */
proc mapimport out=hydro
  datafile="u:\nc_hydrography\hydromaj_arc.shp";
run;


/* This data is all river basins in North Carolina */

/* Fix some of the basin names and delete missing */
data hydro; set hydro;
  label riv_basin='River Basin';
  riv_basin = propcase(riv_basin);
  if riv_basin='' then delete;
run;

/* Limit to only Cape Fear River Basin */
data hydro1; set hydro;
  where riv_basin='Cape Fear';
run;

/* Convert state plane coordinates to WGS 84 coorinates.   */
/* Use http://spatialreference.org using information from  */
/* the .prj file and from sashelp.proj4def.                */
proc gproject data=hydro1 out=hydro1
from="+proj=lcc +lat_1=34.33333333333334 +lat_2=36.16666666666666 
      +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +ellps=GRS80 
      +datum=NAD83 +units=m +no_defs"
    to="EPSG:4326";
    id objectid;
run;

/* Add missing y value for series segments */
proc sort data=hydro1 out=hydro1; by objectid; run;
data hydro1; set hydro1;
  by objectid;
  if last.objectid then do;
    output;
    y=.; output;
  end;
  else output;
  run;
run;

title1 'SGMAP with SERIES Plot, Missing Y Values';
title2 'Cape Fear River and Tributaries';
proc sgmap plotdata=hydro1;
  openstreetmap;
  series x=x y=y / legendlabel='Cape Fear Basin';
run;

/* Hurricane Florence data */
data florence;
  infile datalines dlm=',';
  length lat 8. long 8. wind 8. pres 8. type $32 date $32;
  input lat long wind pres type date;
  datalines;
33.6,76.0,105,955,Hurricane,Sept 13
33.7,76.2,100,955,Hurricane,
33.9,76.4,100,955,Hurricane,
34.0,76.8,90,956,Hurricane,
34.1,77.2,90,954,Hurricane,Sept 14
34.2,77.4,90,958,Hurricane,
34.1,77.9,90,958,Hurricane,
34.0,78.0,80,958,Hurricane,
34.0,78.4,75,968,Hurricane,       
34.0,78.6,70,972,Tropical Storm,
33.9,78.8,70,975,Tropical Storm,
33.8,79.1,65,980,Tropical Storm,
33.7,79.3,60,984,Tropical Storm,Sept 15
33.6,79.5,50,986,Tropical Storm,
33.6,79.5,50,989,Tropical Storm,
33.6,79.6,45,995,Tropical Storm,
33.6,79.8,45,997,Tropical Storm,
33.6,79.9,45,997,Tropical Storm,
33.6,80.1,45,997,Tropical Storm,
33.7,80.5,40,998,Tropical Storm,
33.7,80.8,40,999,Tropical Storm,Sept 16
33.8,81.4,35,1000,Tropical Depression,
34.0,81.8,35,1002,Tropical Depression,
34.6,82.2,35,1006,Tropical Depression,
35.5,82.1,30,1007,Tropical Depression,
36.9,82.2,30,1008,Tropical Depression,Sept 17
;

/* Make the longitude negative  */
data florence; set florence;
  long = -long;
run;

/* Concatenate data sets */
data both; set hydro1 florence; run;

title1 'SGMAP with SERIES and BUBBLE Plots';
title2 'Bubble Size Represents Wind Speed';
proc sgmap plotdata=both;
  openstreetmap;
  /* series plot for rivers */
  series x=x y=y;
  /* plots for the tropical storm */
  series x=long y=lat / smoothconnect;
  bubble x=long y=lat size=wind / transparency=.6
    group=type nomissinggroup
    datalabel=date datalabelpos=bottomright
	datalabelattrs=(color=cxf42941 size=10)
    name="bubble";
  /* add a title to the legend */
  keylegend "bubble" / title="Storm Classification";
run;

quit;
ods _all_ close;