Combine Raster temperature with Mexico municipalities Shapefile

library(readstata13)
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

library(ncdf4)
library(raster)

Loading required package: sp

Attaching package: 'raster'

The following object is masked from 'package:dplyr':

    select

library(sp)
library(rgeos)

rgeos version: 0.3-26, (SVN revision 560)
 GEOS runtime version: 3.5.1-CAPI-1.9.1 r4246 
 Linking to sp version: 1.2-4 
 Polygon checking: TRUE 

library(rgdal)

rgdal: version: 1.2-15, (SVN revision 691)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
 Path to GDAL shared files: /usr/share/gdal/2.2
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.2-4 

library(dismo)

Read ShapeFile and Raster File.

muni.shp<-readOGR(dsn="gis/shapefiles/conjunto_de_datos/",layer = "areas_geoestadisticas_municipales")

OGR data source with driver: ESRI Shapefile 
Source: "gis/shapefiles/conjunto_de_datos/", layer: "areas_geoestadisticas_municipales"
with 2458 features
It has 3 fields

RASTER_FILE_NAME <- 'gis/noaa_temperatures/air.mon.mean.nc'

DATES <- c(836:838)

r.stack<-stack(RASTER_FILE_NAME,bands=DATES)
r.stack<- rotate(r.stack)

data.frame( rasterToPoints( r.stack ) )

mexico_limits <- extent(-118,-82,14,34)
r.stack <- crop(r.stack, mexico_limits)
print (r.stack)

class       : RasterBrick 
dimensions  : 40, 72, 2880, 3  (nrow, ncol, ncell, nlayers)
resolution  : 0.5, 0.5  (x, y)
extent      : -118, -82, 14, 34  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : X2017.08.01.01.06.32, X2017.09.01.01.06.32, X2017.10.01.01.06.32 
min values  :               288.13,               286.52,               281.90 
max values  :               314.63,               311.87,               312.26 

#df.r<- as.data.frame(rc,xy=TRUE,optional = TRUE)
plot(r.stack)

Project both Shapefile and Raster on WGS84 and check that they match on the map.

# tranform both data set to longlat
#myshp<- subset(myshp,CVE_ENT=='01')

rgeo <- projectRaster(r.stack, crs='+proj=longlat +datum=WGS84')
pgeo <- spTransform(muni.shp, CRS('+proj=longlat +datum=WGS84'))

# plot on top of Google map
e <- union(extent(rgeo), extent(pgeo))
print(e)

class       : Extent 
xmin        : -120 
xmax        : -80 
ymin        : 12 
ymax        : 36 

g <- gmap(e,lonlat=TRUE)
plot(g, inter=TRUE)

plot(extent(rgeo), add=TRUE, col='red', lwd=2)
plot(pgeo, add=TRUE, col='blue', lwd=2)

kelvins_to_celsius= function(x){
  return (x-273.15)
}

#extraction <-extract(rgeo,pgeo, weights = TRUE, normalizeWeights=TRUE, fun=mean,df=TRUE)
extraction <-extract(rgeo,pgeo,fun=mean,df=TRUE)

extraction <- extraction  %>% mutate_at(vars(-ID),funs(kelvins_to_celsius))


pgeo@data <- data.frame(pgeo@data, means=extraction)

spplot(pgeo, "means.X2017.08.01.01.06.32",lwd=.1,title ='Mean temp (Celsius) for August 2014')

Read and add tariffa to shapefile data :

pgeo@data <- pgeo@data %>% mutate(id = paste(CVE_ENT,CVE_MUN,sep=''))

df.tariffa<-  read.dta13('gis/tariff/Tariffs.dta')
summary(df.tariffa)

     date2           munid           tarifa            year     
 Min.   :600.0   Min.   : 1001   DAC    :167899   Min.   :2010  
 1st Qu.:620.0   1st Qu.:13072   01     :147096   1st Qu.:2011  
 Median :642.0   Median :20119   1B     : 59076   Median :2013  
 Mean   :641.8   Mean   :19113   1A     : 43276   Mean   :2013  
 3rd Qu.:663.0   3rd Qu.:25013   1C     : 34690   3rd Qu.:2015  
 Max.   :683.0   Max.   :99999   1D     : 12743   Max.   :2016  
                                 (Other): 11808                 
     mun                 id                month       
 Length:476588      Length:476588      Min.   : 1.000  
 Class :character   Class :character   1st Qu.: 4.000  
 Mode  :character   Mode  :character   Median : 7.000  
                                       Mean   : 6.512  
                                       3rd Qu.:10.000  
                                       Max.   :12.000  
                                                       

pgeo@data<- left_join(pgeo@data,df.tariffa, by='id')

df.temp.tariffa.munid<-as.data.frame(pgeo@data)
head(df.temp.tariffa.munid)