In [5]:
setwd("/Users/erinlarson/Dropbox/Documents/GIT/DDESSS_Stream_Dream_Team")

In [7]:
TPvsSO = read.csv("Lat_Long_P_LandUse_StreamOrder.csv")
t(TPvsSO[1,])


1
ActivityStartDate2/7/01
mo2
day7
date2001-02-07T00:00:00Z
MonitoringLocationIdentifierUSGS-09498400
HydrologicConditionStable, normal stage
HydrologicEventRoutine sample
ResultDetectionConditionTextNot Detected
CharacteristicNamePhosphorus
ResultSampleFractionTextTotal
ResultMeasureValueNA
ResultMeasure20
ResultMeasure_MeasureUnitCode
ResultStatusIdentifierAccepted
HUCEightDigitCode15060103
DrainageAreaMeasure_MeasureValue195
DrainageAreaMeasure_MeasureUnitCodesq mi
HUC815060103
PFOR273.918
PWETL20.01
PURB20.132
PAGT20.076
OrganizationIdentifierUSGS-AZ
OrganizationFormalNameUSGS Arizona Water Science Center
MonitoringLocationIdentifier.1USGS-09498400
MonitoringLocationNamePINAL CREEK AT INSPIRATION DAM, NR GLOBE, AZ.
MonitoringLocationTypeNameStream
MonitoringLocationDescriptionTextNA
HUCEightDigitCode.115060103
DrainageAreaMeasure_MeasureValue.1195
DrainageAreaMeasure_MeasureUnitCode.1sq mi
ContributingDrainageAreaMeasure_MeasureValueNA
ContributingDrainageAreaMeasure_MeasureUnitCode
LatitudeMeasure33.57311
LongitudeMeasure-110.9012
SourceMapScaleNumericNA
HorizontalAccuracyMeasure.MeasureValue10
HorizontalAccuracyMeasure.MeasureUnitCodeseconds
HorizontalCollectionMethodNameInterpolated from MAP.
HorizontalCoordinateReferenceSystemDatumNameNAD83
VerticalMeasure.MeasureValue2740
VerticalMeasure.MeasureUnitCodefeet
VerticalAccuracyMeasure.MeasureValue20
VerticalAccuracyMeasure.MeasureUnitCodefeet
VerticalCollectionMethodNameInterpolated from topographic map.
VerticalCoordinateReferenceSystemDatumNameNGVD29
CountryCodeUS
StateCode4
CountyCode7
AquiferNameNA
FormationTypeTextNA
AquiferTypeNameNA
ConstructionDateTextNA
WellDepthMeasure.MeasureValueNA
WellDepthMeasure.MeasureUnitCodeNA
WellHoleDepthMeasure.MeasureValueNA
WellHoleDepthMeasure.MeasureUnitCodeNA
ProviderNameNWIS
HUC8.115060103
STRAHLER22.4

In [8]:
TPvsSO2<-subset(TPvsSO, ResultMeasure2<4.0, select=c(ResultMeasure2, STRAHLER2, DrainageAreaMeasure_MeasureValue, HUCEightDigitCode, date, ResultSampleFractionText, PFOR2, PWETL2,PURB2, PAGT2, LatitudeMeasure, LongitudeMeasure, MonitoringLocationIdentifier))
head(TPvsSO2)
write.csv(TPvsSO2, "WorkingDataSet.csv")


ResultMeasure2STRAHLER2DrainageAreaMeasure_MeasureValueHUCEightDigitCodedateResultSampleFractionTextPFOR2PWETL2PURB2PAGT2LatitudeMeasureLongitudeMeasureMonitoringLocationIdentifier
0 2.4 195 15060103 2001-02-07T00:00:00ZTotal 73.918 0.01 0.132 0.076 33.57311 -110.9012 USGS-09498400
0 2.4 195 15060103 2001-04-19T00:00:00ZTotal 73.918 0.01 0.132 0.076 33.57311 -110.9012 USGS-09498400
0 2.4 195 15060103 2001-09-04T00:00:00ZTotal 73.918 0.01 0.132 0.076 33.57311 -110.9012 USGS-09498400
0 2.4 195 15060103 2001-12-20T00:00:00ZTotal 73.918 0.01 0.132 0.076 33.57311 -110.9012 USGS-09498400
0 2.4 195 15060103 2002-06-11T00:00:00ZTotal 73.918 0.01 0.132 0.076 33.57311 -110.9012 USGS-09498400
0 2.4 195 15060103 2002-09-06T00:00:00ZTotal 73.918 0.01 0.132 0.076 33.57311 -110.9012 USGS-09498400

In [9]:
summary(TPvsSO)


 ActivityStartDate       mo              day                          date     
 7/17/01:  18      Min.   : 1.000   Min.   : 1.00   2001-07-17T00:00:00Z:  18  
 4/17/07:  17      1st Qu.: 4.000   1st Qu.:10.00   2007-04-17T00:00:00Z:  17  
 6/17/08:  17      Median : 7.000   Median :17.00   2008-06-17T00:00:00Z:  17  
 3/29/00:  15      Mean   : 6.534   Mean   :16.54   2000-03-29T00:00:00Z:  15  
 6/13/01:  15      3rd Qu.: 9.000   3rd Qu.:22.00   2000-07-19T00:00:00Z:  15  
 7/19/00:  15      Max.   :12.000   Max.   :31.00   2001-06-13T00:00:00Z:  15  
 (Other):4955                                       (Other)             :4955  
 MonitoringLocationIdentifier           HydrologicCondition
 USGS-09041400: 250           Stable, normal stage:5052    
 USGS-09041090: 230                                        
 USGS-09025010: 228                                        
 USGS-09027100: 228                                        
 USGS-09023750: 227                                        
 USGS-09013500: 199                                        
 (Other)      :3690                                        
       HydrologicEvent ResultDetectionConditionText  CharacteristicName
 Routine sample:5052               :4575            Phosphorus:5052    
                       Not Detected: 477                               
                                                                       
                                                                       
                                                                       
                                                                       
                                                                       
 ResultSampleFractionText ResultMeasureValue ResultMeasure2    
 Dissolved:1349           Min.   : 0.0020    Min.   : 0.00000  
 Total    :3703           1st Qu.: 0.0110    1st Qu.: 0.00800  
                          Median : 0.0250    Median : 0.02100  
                          Mean   : 0.0812    Mean   : 0.07354  
                          3rd Qu.: 0.0580    3rd Qu.: 0.05100  
                          Max.   :24.0000    Max.   :24.00000  
                          NA's   :477                          
 ResultMeasure_MeasureUnitCode ResultStatusIdentifier HUCEightDigitCode 
          : 477                Accepted   :4355       Min.   :14010001  
 mg/l as P:4575                Historical : 277       1st Qu.:14010001  
                               Preliminary: 420       Median :14010001  
                                                      Mean   :14155881  
                                                      3rd Qu.:14040106  
                                                      Max.   :15060203  
                                                                        
 DrainageAreaMeasure_MeasureValue DrainageAreaMeasure_MeasureUnitCode
 Min.   :   1.67                  sq mi:5052                         
 1st Qu.:  46.90                                                     
 Median : 167.00                                                     
 Mean   :1007.27                                                     
 3rd Qu.: 567.00                                                     
 Max.   :7986.00                                                     
                                                                     
      HUC8              PFOR2           PWETL2            PURB2        
 Min.   :14010001   Min.   : 0.00   Min.   :0.00000   Min.   :0.00000  
 1st Qu.:14010001   1st Qu.:46.08   1st Qu.:0.00000   1st Qu.:0.00000  
 Median :14010001   Median :46.08   Median :0.00000   Median :0.00000  
 Mean   :14155881   Mean   :53.84   Mean   :0.02283   Mean   :0.01712  
 3rd Qu.:14040106   3rd Qu.:72.01   3rd Qu.:0.00000   3rd Qu.:0.00000  
 Max.   :15060203   Max.   :98.94   Max.   :4.00000   Max.   :0.50667  
                                                                       
     PAGT2        OrganizationIdentifier
 Min.   :0.0000   USGS-AZ: 451          
 1st Qu.:0.0000   USGS-CO:4288          
 Median :0.2500   USGS-NM:  90          
 Mean   :0.2702   USGS-NV: 132          
 3rd Qu.:0.2500   USGS-UT:  11          
 Max.   :3.0500   USGS-WY:  80          
                                        
                          OrganizationFormalName MonitoringLocationIdentifier.1
 USGS Arizona Water Science Center   : 451       USGS-09041400: 250            
 USGS Colorado Water Science Center  :4288       USGS-09041090: 230            
 USGS Nevada Water Science Center    : 132       USGS-09025010: 228            
 USGS New Mexico Water Science Center:  90       USGS-09027100: 228            
 USGS Utah Water Science Center      :  11       USGS-09023750: 227            
 USGS Wyoming Water Science Center   :  80       USGS-09013500: 199            
                                                 (Other)      :3690            
                                        MonitoringLocationName
 MUDDY CRK BLW WOLFORD MTN RESER. NR KREMMLING, CO : 250      
 MUDDY CREEK ABOVE ANTELOPE CREEK NR. KREMMLING, CO: 230      
 FRASER RIVER AT TABERNASH, CO                     : 228      
 FRASER RIVER BLW VASQUEZ CREEK AT WINTER PARK CO. : 228      
 FRASER RIVER BELOW BUCK CREEK AT WINTER PARK, CO  : 227      
 EAST INLET NEAR GRAND LAKE, CO.                   : 199      
 (Other)                                           :3690      
 MonitoringLocationTypeName MonitoringLocationDescriptionText
 Stream:5052                Mode:logical                     
                            NA's:5052                        
                                                             
                                                             
                                                             
                                                             
                                                             
 HUCEightDigitCode.1 DrainageAreaMeasure_MeasureValue.1
 Min.   :14010001    Min.   :   1.67                   
 1st Qu.:14010001    1st Qu.:  46.90                   
 Median :14010001    Median : 167.00                   
 Mean   :14155881    Mean   :1007.27                   
 3rd Qu.:14040106    3rd Qu.: 567.00                   
 Max.   :15060203    Max.   :7986.00                   
                                                       
 DrainageAreaMeasure_MeasureUnitCode.1
 sq mi:5052                           
                                      
                                      
                                      
                                      
                                      
                                      
 ContributingDrainageAreaMeasure_MeasureValue
 Min.   :  10.5                              
 1st Qu.:  59.0                              
 Median :  59.0                              
 Mean   :1069.2                              
 3rd Qu.:2829.0                              
 Max.   :5494.0                              
 NA's   :4427                                
 ContributingDrainageAreaMeasure_MeasureUnitCode LatitudeMeasure
      :4427                                      Min.   :32.65  
 sq mi: 625                                      1st Qu.:38.88  
                                                 Median :39.99  
                                                 Mean   :39.19  
                                                 3rd Qu.:40.14  
                                                 Max.   :43.03  
                                                                
 LongitudeMeasure SourceMapScaleNumeric HorizontalAccuracyMeasure.MeasureValue
 Min.   :-113.9   Min.   :24000         0.1    :  51                          
 1st Qu.:-107.6   1st Qu.:24000         0.5    : 496                          
 Median :-106.4   Median :24000         1      :2579                          
 Mean   :-107.2   Mean   :24025         10     : 408                          
 3rd Qu.:-105.8   3rd Qu.:24000         5      :1472                          
 Max.   :-105.7   Max.   :25000         Unknown:  46                          
                  NA's   :1861                                                
 HorizontalAccuracyMeasure.MeasureUnitCode
 minutes: 365                             
 seconds:4641                             
 Unknown:  46                             
                                          
                                          
                                          
                                          
                                              HorizontalCollectionMethodName
 Interpolated from Digital MAP.                              : 698          
 Interpolated from MAP.                                      :4248          
 Mapping grade GPS unit (handheld accuracy range 12 to 40 ft): 106          
                                                                            
                                                                            
                                                                            
                                                                            
 HorizontalCoordinateReferenceSystemDatumName VerticalMeasure.MeasureValue
 NAD83:5052                                   Min.   :1570                
                                              1st Qu.:6100                
                                              Median :8006                
                                              Mean   :7051                
                                              3rd Qu.:8380                
                                              Max.   :9960                
                                              NA's   :61                  
 VerticalMeasure.MeasureUnitCode VerticalAccuracyMeasure.MeasureValue
     :  61                       Min.   : 1.00                       
 feet:4991                       1st Qu.:10.00                       
                                 Median :15.00                       
                                 Mean   :13.98                       
                                 3rd Qu.:20.00                       
                                 Max.   :50.00                       
                                 NA's   :61                          
 VerticalAccuracyMeasure.MeasureUnitCode
     :  61                              
 feet:4991                              
                                        
                                        
                                        
                                        
                                        
                            VerticalCollectionMethodName
                                          :  61         
 Interpolated from Digital Elevation Model: 188         
 Interpolated from topographic map.       :3256         
 Level or other surveyed method.          : 386         
 Unknown.                                 :1161         
                                                        
                                                        
 VerticalCoordinateReferenceSystemDatumName CountryCode   StateCode     
       :  61                                US:5052     Min.   : 4.000  
 NAVD88: 280                                            1st Qu.: 8.000  
 NGVD29:4711                                            Median : 8.000  
                                                        Mean   : 8.906  
                                                        3rd Qu.: 8.000  
                                                        Max.   :56.000  
                                                                        
   CountyCode     AquiferName    FormationTypeText AquiferTypeName
 Min.   :  3.00   Mode:logical   Mode:logical      Mode:logical   
 1st Qu.: 49.00   NA's:5052      NA's:5052         NA's:5052      
 Median : 49.00                                                   
 Mean   : 50.32                                                   
 3rd Qu.: 49.00                                                   
 Max.   :107.00                                                   
                                                                  
 ConstructionDateText WellDepthMeasure.MeasureValue
 Mode:logical         Mode:logical                 
 NA's:5052            NA's:5052                    
                                                   
                                                   
                                                   
                                                   
                                                   
 WellDepthMeasure.MeasureUnitCode WellHoleDepthMeasure.MeasureValue
 Mode:logical                     Mode:logical                     
 NA's:5052                        NA's:5052                        
                                                                   
                                                                   
                                                                   
                                                                   
                                                                   
 WellHoleDepthMeasure.MeasureUnitCode ProviderName     HUC8.1        
 Mode:logical                         NWIS:5052    Min.   :14010001  
 NA's:5052                                         1st Qu.:14010001  
                                                   Median :14010001  
                                                   Mean   :14155881  
                                                   3rd Qu.:14040106  
                                                   Max.   :15060203  
                                                                     
   STRAHLER2    
 Min.   :1.000  
 1st Qu.:2.000  
 Median :2.000  
 Mean   :2.272  
 3rd Qu.:2.000  
 Max.   :7.000  
                

In [12]:
library(ggplot2)
library(lattice)
library(RColorBrewer)


Warning message:
“package ‘ggplot2’ was built under R version 3.3.2”Warning message:
“package ‘lattice’ was built under R version 3.3.2”

In [13]:
install.packages("ggplot2")


The downloaded binary packages are in
	/var/folders/7d/vkx258d969q4tfdk2vgfff0r0000gn/T//RtmpKDWZqY/downloaded_packages

In [14]:
install.packages("RColorBrewer")


The downloaded binary packages are in
	/var/folders/7d/vkx258d969q4tfdk2vgfff0r0000gn/T//RtmpKDWZqY/downloaded_packages

In [15]:
install.packages("lattice")


The downloaded binary packages are in
	/var/folders/7d/vkx258d969q4tfdk2vgfff0r0000gn/T//RtmpKDWZqY/downloaded_packages

In [16]:
display.brewer.all()


TDP Time series with dataset with outliers removed


In [35]:
TDP.HUC.nooutlier.timeplot<-xyplot(ResultMeasure2~date|HUCEightDigitCode, 
                     data=subset(TPvsSO2, ResultSampleFractionText=="Dissolved"), 
                     type="a", xlab="Date", ylab="Total Dissolved Phosphorus mg/L")
TDP.HUC.nooutlier.timeplot


Total P Time series with dataset with outliers removed


In [36]:
TotalP.HUC.nooutlier.timeplot<-xyplot(ResultMeasure2~date|HUCEightDigitCode, 
                     data=subset(TPvsSO2, ResultSampleFractionText=="Total"), 
                     type="a", xlab="Date", ylab="Total Phosphorus mg/L")
TotalP.HUC.nooutlier.timeplot


Total P Time Series by Site


In [34]:
TotalP.Site.nooutlier.timeplot<-xyplot(ResultMeasure2~date|MonitoringLocationIdentifier, 
                     data=subset(TPvsSO2, ResultSampleFractionText=="Total"), 
                     type="a", xlab="Date", ylab="Total Phosphorus mg/L")

TotalP.Site.nooutlier.timeplot
trellis.device(device="png", filename="TotalPTimeSeriesbySiteNoOutliers.png")


TDP Time Series by Site


In [33]:
TDP.Site.nooutlier.timeplot<-xyplot(ResultMeasure2~date|MonitoringLocationIdentifier, 
                     data=subset(TPvsSO2, ResultSampleFractionText=="Dissolved"), 
                     type="a", xlab="Date", ylab="Total Dissolved Phosphorus mg/L")

TDP.Site.nooutlier.timeplot


TP by Stream Order


In [42]:
TP.SO.nooutlier<-plot(ResultMeasure2~STRAHLER2, 
                     data=subset(TPvsSO2, ResultSampleFractionText=="Total"), 
                     type="p", xlab="Stream Order", ylab="Total Phosphorus mg/L")

TP.SO.nooutlier


NULL

TDP by Stream Order


In [43]:
TDP.SO.nooutlier<-plot(ResultMeasure2~STRAHLER2, 
                     data=subset(TPvsSO2, ResultSampleFractionText=="Dissolved"), 
                     type="p", xlab="Stream Order", ylab="Total Dissolved Phosphorus mg/L")

TDP.SO.nooutlier


NULL

Standard deviation by site


In [44]:
library("dplyr")


Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


In [53]:
TotalP <-subset(TPvsSO2, ResultSampleFractionText=="Total")

TP.sd <- group_by(TotalP, MonitoringLocationIdentifier) %>% 
    summarise(sd_Total_P = sd(ResultMeasure2, na.rm = TRUE),
             STRAHLER2 = mean(STRAHLER2, na.rm = TRUE))

In [54]:
head(TP.sd)


MonitoringLocationIdentifiersd_Total_PSTRAHLER2
USGS-09010500 NA 2
USGS-090110000.029114215 2
USGS-090124000.000000000 2
USGS-090125000.011261010 2
USGS-090135000.027024594 2
USGS-090140000.006311886 2

In [55]:
TDP <-subset(TPvsSO2, ResultSampleFractionText=="Dissolved")


TDP.sd <- group_by(TDP, MonitoringLocationIdentifier) %>% 
    summarise(sd_TDP = sd(ResultMeasure2, na.rm = TRUE),
             STRAHLER2 = mean(STRAHLER2, na.rm = TRUE))

In [56]:
head(TDP.sd)


MonitoringLocationIdentifiersd_TDPSTRAHLER2
USGS-090110000.006136319 2
USGS-090125000.005185038 2
USGS-090135000.003677473 2
USGS-09014000 NA 2
USGS-090165000.002486729 2
USGS-090180000.023735346 2

Standard deviation by stream order


In [87]:
variationTP <- plot(log(sd_Total_P+1)~STRAHLER2, data=TP.sd, 
                     type="p", xlab="Stream Order", ylab="Log(TP Standard Deviation by site+1)")

variationTP


NULL

In [86]:
variationTDP <- plot(log(sd_TDP+1)~STRAHLER2, data=TDP.sd, 
                     type="p", xlab="Stream Order", ylab="log(TDP Standard Deviation by site+1)")

variationTDP


NULL

Histogram of sd


In [63]:
hist(TDP.sd$sd_TDP)



In [64]:
hist(log(TDP.sd$sd_TDP))



In [66]:
hist(log(TP.sd$sd_Total_P))


Regression Stream Order vs SD in TP


In [80]:
TP.sd$log_TP_sd<-log(TP.sd$sd_Total_P+1)
TP.sd.no.na<-na.omit(TP.sd)

summary(TP.sd.no.na)


 MonitoringLocationIdentifier   sd_Total_P        STRAHLER2   
 USGS-09011000: 1             Min.   :0.00000   Min.   :1.00  
 USGS-09012400: 1             1st Qu.:0.01535   1st Qu.:2.00  
 USGS-09012500: 1             Median :0.03375   Median :2.00  
 USGS-09013500: 1             Mean   :0.08105   Mean   :2.47  
 USGS-09014000: 1             3rd Qu.:0.07862   3rd Qu.:3.00  
 USGS-09015000: 1             Max.   :0.52297   Max.   :7.00  
 (Other)      :74                                             
   log_TP_sd      
 Min.   :0.00000  
 1st Qu.:0.01523  
 Median :0.03320  
 Mean   :0.07277  
 3rd Qu.:0.07568  
 Max.   :0.42066  
                  

In [81]:
TP.sd.lm <- lm(log_TP_sd~STRAHLER2, data=TP.sd.no.na)
summary(TP.sd.lm)


Call:
lm(formula = log_TP_sd ~ STRAHLER2, data = TP.sd.no.na)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.07371 -0.05879 -0.03803  0.00471  0.34777 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)  0.077073   0.028351   2.719  0.00808 **
STRAHLER2   -0.001743   0.010558  -0.165  0.86928   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0994 on 78 degrees of freedom
Multiple R-squared:  0.0003494,	Adjusted R-squared:  -0.01247 
F-statistic: 0.02726 on 1 and 78 DF,  p-value: 0.8693

Regression Stream Order vs SD in TDP


In [82]:
TDP.sd$log_TDP_sd<-log(TDP.sd$sd_TDP+1)
TDP.sd.no.na<-na.omit(TDP.sd)

summary(TDP.sd.no.na)


 MonitoringLocationIdentifier     sd_TDP           STRAHLER2    
 USGS-09011000: 1             Min.   :0.000000   Min.   :1.000  
 USGS-09012500: 1             1st Qu.:0.003604   1st Qu.:2.000  
 USGS-09013500: 1             Median :0.007599   Median :2.000  
 USGS-09016500: 1             Mean   :0.025923   Mean   :2.313  
 USGS-09018000: 1             3rd Qu.:0.022607   3rd Qu.:3.000  
 USGS-09019000: 1             Max.   :0.443058   Max.   :7.000  
 (Other)      :41                                               
   log_TDP_sd      
 Min.   :0.000000  
 1st Qu.:0.003598  
 Median :0.007571  
 Mean   :0.023870  
 3rd Qu.:0.022356  
 Max.   :0.366764  
                   

In [83]:
TDP.sd.lm <- lm(log_TDP_sd~STRAHLER2, data=TDP.sd.no.na)
summary(TDP.sd.lm)


Call:
lm(formula = log_TDP_sd ~ STRAHLER2, data = TDP.sd.no.na)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.02375 -0.02015 -0.01662 -0.00139  0.34302 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0229737  0.0215294   1.067    0.292
STRAHLER2   0.0003875  0.0085808   0.045    0.964

Residual standard error: 0.05712 on 45 degrees of freedom
Multiple R-squared:  4.531e-05,	Adjusted R-squared:  -0.02218 
F-statistic: 0.002039 on 1 and 45 DF,  p-value: 0.9642