CalidadDelAireSonora.R

#bibliotecas a usar 
library(ggplot2)
library(readxl)
library(RSAGA)

## Loading required package: gstat

## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo

## Loading required package: shapefiles

## Loading required package: foreign

## 
## Attaching package: 'shapefiles'

## The following objects are masked from 'package:foreign':
## 
##     read.dbf, write.dbf

## Loading required package: plyr

library(gstat)
library(rgdal)

## Loading required package: sp

## rgdal: version: 1.4-8, (SVN revision 845)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/mirol/OneDrive/Documentos/R/win-library/3.6/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/mirol/OneDrive/Documentos/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.3-2

library(sf)

## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

library(fdth)

## 
## Attaching package: 'fdth'

## The following objects are masked from 'package:stats':
## 
##     sd, var

library(plotly)

## 
## Attaching package: 'plotly'

## The following objects are masked from 'package:plyr':
## 
##     arrange, mutate, rename, summarise

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(DT)
library(extrafont)

## Registering fonts with R

#Script para el procesamiento y an攼㸱lisis de datos de calidad del aire

#Los datos usados fueron obtenidos de la estaci昼㸳n ubicada 
#en el Instituto Tecnol昼㸳gico de Nogales (ITN):
Nogales <- st_read("ITN.shp")

## Reading layer `ITN' from data source `C:\Users\mirol\OneDrive\Documentos\ProbICIAM\Calidad del aire\ITN.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 504515.8 ymin: 3461860 xmax: 504515.8 ymax: 3461860
## epsg (SRID):    32612
## proj4string:    +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs

mapview::mapview(Nogales ,labels=F)

#Importar datos de PM 10
nogales <- read_excel("nogales.xlsx")
View(nogales)

#Tabla de datos de la estaci昼㸳n de Nogales
datatable(nogales)

#Gr攼㸱fica de datos disponibles en serie de tiempos
nog <- ggplot(nogales, aes(as.Date(nogales$Fecha), nogales$Valor))+ geom_point() + geom_line(colour = 'blue')+xlab("A昼㸱o")+ ylab("Part攼㹤culas PM10")+ggtitle("Calidad del aire estaci昼㸳n Nogales")
ggplotly(nog)

#Gr攼㸱fico de caja y bigote de los datos PM 10
P10 <- ggplot(nogales, aes(x = Fecha, y = Valor)) +
  geom_boxplot(fill = '#99FFCC') + theme(text = element_text(size = 10)) + ggtitle('Calidad de aire Nogales') + theme(plot.title = element_text(size = rel(2), vjust = 2, face = 'bold'))
P10

## Warning: Continuous x aesthetic -- did you forget aes(group=...)?

#Resumen estad攼㹤stico de los datos PM10
summary(nogales)

##      Fecha                         Valor       
##  Min.   :2018-12-05 00:00:00   Min.   :  0.00  
##  1st Qu.:2019-03-04 00:00:00   1st Qu.: 17.00  
##  Median :2019-06-03 00:00:00   Median : 29.00  
##  Mean   :2019-06-03 14:43:57   Mean   : 42.19  
##  3rd Qu.:2019-09-01 00:00:00   3rd Qu.: 52.00  
##  Max.   :2019-12-05 00:00:00   Max.   :445.00

#histograma de los datos PM10
histogramaN <- ggplot(nogales, aes(x=Valor)) + 
  geom_histogram() 
ggplotly(histogramaN)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#Tabla de distribuciones de frecuencias PM10
dist <- fdt(nogales,breaks="Sturges") 
dist

## Valor 
##  Class limits    f   rf rf(%)   cf  cf(%)
##        [0,30) 4367 0.51 50.87 4367  50.87
##     [30,59.9) 2543 0.30 29.62 6910  80.49
##   [59.9,89.9)  799 0.09  9.31 7709  89.80
##    [89.9,120)  372 0.04  4.33 8081  94.13
##     [120,150)  209 0.02  2.43 8290  96.56
##     [150,180)  121 0.01  1.41 8411  97.97
##     [180,210)   64 0.01  0.75 8475  98.72
##     [210,240)   51 0.01  0.59 8526  99.31
##     [240,270)   19 0.00  0.22 8545  99.53
##     [270,300)   14 0.00  0.16 8559  99.70
##     [300,330)   11 0.00  0.13 8570  99.83
##     [330,360)    8 0.00  0.09 8578  99.92
##     [360,390)    4 0.00  0.05 8582  99.97
##     [390,419)    0 0.00  0.00 8582  99.97
##     [419,449)    3 0.00  0.03 8585 100.00

#############################################################################

#Script para el procesamiento y an攼㸱lisis de datos de calidad del aire

HILLO <- st_read("HMO.shp")

## Reading layer `HMO' from data source `C:\Users\mirol\OneDrive\Documentos\ProbICIAM\Calidad del aire\HMO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 11 fields
## geometry type:  POINT
## dimension:      XYZ
## bbox:           xmin: 502871.4 ymin: 3217010 xmax: 502871.4 ymax: 3217010
## epsg (SRID):    32612
## proj4string:    +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs

mapview::mapview(HILLO, labels=F)

#Importar datos de PM 10
hermosillo <- read_excel("hermosillo.xlsx")
View(hermosillo)

#Tabla de datos de la estaci昼㸳n de Nogales
datatable(hermosillo)

#Gr攼㸱fica de datos disponibles en serie de tiempos
hillo <- ggplot(hermosillo, aes(as.Date(hermosillo$Fecha), hermosillo$PM10))+ geom_point() + geom_line(colour = 'tomato')+ xlab("A昼㸱o")+ ylab("Part攼㹤culas PM10")+ggtitle("Calidad del aire estaci昼㸳n Hermosillo")
ggplotly(hillo)

#Gr攼㸱fico de caja y bigote de los datos PM 10
PM10 <- ggplot(hermosillo, aes(x = Fecha, y = PM10)) +
  geom_boxplot(fill = '#FFCC99') + theme(text = element_text(size = 10)) + ggtitle('Calidad de aire Hermosillo') + theme(plot.title = element_text(size = rel(2), vjust = 2, face = 'bold'))
PM10

## Warning: Continuous x aesthetic -- did you forget aes(group=...)?

#Resumen estad攼㹤stico de los datos PM10
summary(hermosillo)

##      Fecha                          PM10       
##  Min.   :2019-01-01 00:00:00   Min.   :  4.07  
##  1st Qu.:2019-03-19 00:00:00   1st Qu.: 18.69  
##  Median :2019-05-31 12:00:00   Median : 27.40  
##  Mean   :2019-06-02 00:09:50   Mean   : 31.77  
##  3rd Qu.:2019-08-16 00:00:00   3rd Qu.: 39.43  
##  Max.   :2019-11-03 00:00:00   Max.   :188.22

#histograma de los datos PM10
histogramaH <- ggplot(hermosillo, aes(x=PM10)) + geom_histogram() 
ggplotly(histogramaH)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#Tabla de distribuciones de frecuencias PM10
dist <- fdt(hermosillo,breaks="Sturges") 
dist

## PM10 
##       Class limits    f   rf rf(%)   cf  cf(%)
##   [4.0293,17.3202) 1522 0.22 21.69 1522  21.69
##  [17.3202,30.6111) 2544 0.36 36.25 4066  57.94
##  [30.6111,43.9021) 1588 0.23 22.63 5654  80.56
##   [43.9021,57.193)  709 0.10 10.10 6363  90.67
##   [57.193,70.4839)  339 0.05  4.83 6702  95.50
##  [70.4839,83.7748)  142 0.02  2.02 6844  97.52
##  [83.7748,97.0658)   68 0.01  0.97 6912  98.49
##  [97.0658,110.357)   40 0.01  0.57 6952  99.06
##  [110.357,123.648)   26 0.00  0.37 6978  99.43
##  [123.648,136.939)   22 0.00  0.31 7000  99.74
##  [136.939,150.229)   14 0.00  0.20 7014  99.94
##   [150.229,163.52)    3 0.00  0.04 7017  99.99
##   [163.52,176.811)    0 0.00  0.00 7017  99.99
##  [176.811,190.102)    1 0.00  0.01 7018 100.00

#box plot comparativo
boxplot(nogales$Valor, hermosillo$PM10, horizontal = TRUE, col = c('#99FFCC', '#FFCC99'), names = c('Nogales', 'Hermosillo'), xlab='PM10', main = 'Comparacion de las dos zonas')