In IREA-CNR-MI/sprawl: sprawl: Spatial Processing (in) R: Advanced Workflows Library

knitr::opts_chunk$set(echo = TRUE)

Results of Phenological Mapping - ERMES project

This document show results of rice phenological mapping conducted within the ERMES project for the r params$country study area. Sowing dates were estimated starting from time series of MODIS 250m resolution imagery, using the PhenoRice algorithm.

Raster Maps

Raster data - aggregated on 2x2 Km regular grid

2003-2016 Sowing dates maps

library(raster)
library(xts)
library(sp)
library(gdalUtils)
library(rgdal)
library(data.table)
library(dplyr)
library(plyr)
library(ggplot2)
library(reshape)
library(grid)
library(gridExtra)
library(hash)
library("scales")
library(tools)
library(stringr)
library(knitr)
library(maptools)
library(ireaRscripts)
library(lubridate)
library(RPostgreSQL)
library(RgoogleMaps)
library(ggmap)
# countries = c('IT','ES','GR')
# country = c('IT')
country = params$country
cc = country

Main_Folder = "/home/lb/projects/ermes/datasets/rs_products/Phenology/%cc%/2016/v1.0/Outputs/ERMES_Grid"
Grid_Folder = "/home/lb/projects/ermes/datasets/ERMES_Folder_Structure/%cc%/Regional/%cc%_Reference_Grid/%cc%_ERMES_Regional_Grid.shp"
admin_shape = "/home/lb/projects/ermes/datasets/rs_products/Phenology/Ancillary_Datasets/World_Provinces/provinces_world_laea_ermes.shp"
out_folder = file.path(str_replace_all(Main_Folder,"%cc%",country),'Summaries','pdf_and_graphs')
out_folder_anomalies = file.path(str_replace_all(Main_Folder,"%cc%",country),'Summaries','Anomalies')

vers = '1.0'
selvar = c(1,1,1,1)

start_year = start_sel_years = 2003
end_year = end_sel_years = 2016
laea_crs = CRS("+init=epsg:3035 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0 +units=m +no_defs")
geo_WGS84_crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0,0,0,0 ")

ermes_grid = str_replace_all(Grid_Folder,"%cc%",country)
in_RData_file = file.path(str_replace_all(Main_Folder,"%cc%",country),'RData',paste0('ERMES_Aggregates_',start_year,'_',end_year,'.RData'))
data_in = get(load(file = in_RData_file))
selected = as.character(seq(start_sel_years,end_sel_years, 1))
data_in = droplevels(data_in[data_in$year %in% selected,])  # OKKIO !!!!

# Retieve spatialpointsdataframe of cells for selected country

ERMES_cells_poly = readOGR(dirname(ermes_grid) , basename(file_path_sans_ext(ermes_grid)))
ext_grd = extent(ERMES_cells_poly)
# ext_grd_gg = spTransform(ext_grd,CRSobj = CRS("+init=epsg:4326"))

#  Retieve a polygon map of italian regions

mapit = readOGR(dirname(admin_shape),file_path_sans_ext(basename(admin_shape)), drop_unsupported_fields=T)
mapit = spTransform(mapit, geo_WGS84_crs)
mapit_df = fortify(mapit, region = "name")
mapit_data = data.frame(id = unique(mapit_df$id), value = rnorm(174))

# Add spatial info to the data_in fata frame + add some dummy variables
data_in$group = 1
data_in = join(data_in, ERMES_cells_poly@data, by = 'int_id')
data_in$variable = as.factor(data_in$variable)                 #  convert "variable" column to factor
data_in$year = as.factor(data_in$year)                 #  convert "year" column to factor
is.na(data_in) <- do.call(cbind,lapply(data_in, is.nan))
data_in = subset(data_in, !is.na(mean))                                                                 # remove data with "NA" in the mean column
data_in[,"Date":=as.Date(strptime(paste("2008", data_in$mean), format="%Y %j") ),with=FALSE]
data_in$percol <- cut(100*data_in$perc_rice, breaks = c(0,1,10,20,30,40,50,60,70,80,90,100,110))     # catgorize the rice fc - 10 classes
data_in$variable = factor(data_in$variable,levels(data_in$variable)[c(3,4,1,2)])  # reorder the variables
data_in$rice_area = data_in$perc_rice*2000*2000           # compute retrieved area

data_in = subset(data_in, is_rice == 1)  # Consider only cells of "Rice" muncipalities
data_in$variable = factor(data_in$variable, levels = c( 'MinDoys','SoSDoys', 'MaxDoys', 'MaxVis' ))
levels(data_in$variable) = c( 'DOY of Sowing', 'DOY of Emergence', 'DOY of Heading', 'Value of VI at Heading')



mylabels = function(x){

  labels = lb_doytodate(x, 2003)
  labels = format(labels, "%b-%d")
  labels

}

# Build the plot for the rice fraction : Add points to the map, set colors and variables and set limits ----
ncols = 4

data_mindoy = droplevels(subset(data_in, variable == 'DOY of Sowing' & data_in$year %in% selected ))
 data_mindoy$mean[data_mindoy$Rice_fc<0.05] = NA
# Convert to Geographic projection and create a temporary dataframe containing data to be used for plottuing

newdf = NULL
for (yy in selected) {
  yy
  data_temp = droplevels(subset(data_mindoy, year == yy))
  temprast = rasterFromXYZ(data.frame(x = data_temp$x_LAEA, y = data_temp$y_LAEA, z  =data_temp$mean),crs = laea_crs)
  temprast2 = projectRaster(temprast, crs = geo_WGS84_crs, res = c(0.0254, 0.0254))
  temp_df = data.frame(lat = coordinates(temprast2)[,1], lon = coordinates(temprast2)[,2], z = getValues(temprast2))
  temp_df$cellcode = seq(1:length(temp_df$lat))
  temp_df$year = yy
  newdf = rbind(newdf, temp_df)

}

newdf$year = as.factor(newdf$year)
newdf$cellcode = as.factor(newdf$cellcode)


# Computye tha anomalies on the temporarty dataframe using dplyr ----
  grouped = group_by(newdf, cellcode)  
  newdf_avgs =  mutate(grouped, anomaly = (z - mean(z, na.rm = T)))


# Do the plots of sowing dates (multiple years) ----

if (cc == 'IT'){
  gc <- geocode("pavia, it")
  center <- as.numeric(gc)
  xlims = c(8,10)
  ylims = c(44.7,45.8)
  zoom = 8
}

if (cc == 'ES'){
  gc <- geocode("castellon, es")
  center <- as.numeric(gc)
  xlims = c(-0.5,1)
  ylims = c(39,40.8)
  zoom = 8
}

if (cc == 'GR'){
  gc <- geocode("thesasloniki, gr")
  center <- as.numeric(gc)
  xlims = c(22.4, 23.55)
  ylims = c(40.4,41.2)
  zoom = 9
}

G <- ggmap(get_map(location = center,  zoom = zoom, scale = 2,maptype = 'terrain', source = 'google')) + xlim(xlims) + ylim(ylims)
mapmin = G + geom_tile(data =  newdf,aes(x = lat, y = lon, fill = z), width = 0.0254 ,height = 0.0254)
mapmin <- mapmin + scale_fill_gradientn('Date of Sowing', limits=c(90, 160),colours = RColorBrewer::brewer.pal(5,"RdYlGn"), labels = mylabels, na.value = 'transparent', guide = "legend", breaks = seq(90,160,10))  
if (country == "IT"){ mapmin = mapmin +coord_fixed(xlim = c(8,10))}
mapmin = mapmin +coord_fixed(xlim = xlims)
mapmin = mapmin + facet_wrap(~year)
mapmin = mapmin + theme_bw() + labs(title = "Estimated dates of Sowing", x = "Longitude", y = "Latitude") +
  theme(plot.title = element_text(size = 14, vjust = 1)) + labs(x = "Longitude") +
  theme(axis.text.x  = element_text(size = 8), axis.text.y  = element_text(size = 8))+
  theme(legend.position="right")+theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), 'cm'))+
  theme(legend.background = element_rect(colour = 'black'))+
  theme(legend.justification = c(0,0) ,legend.background = element_rect(colour = 'black'))
print(mapmin)

2016 Sowing dates map

newdf_sub =  droplevels(subset(newdf,year %in% 2016 ))

G <- ggmap(get_map(location = center,  zoom = zoom, scale = 2,maptype = 'terrain', source = 'google')) + xlim(xlims) + ylim(ylims)
mapmin = G + geom_tile(data =  newdf_sub,aes(x = lat, y = lon, fill = z), width = 0.0254 ,height = 0.0254)
mapmin <- mapmin + scale_fill_gradientn('Date of Sowing', limits=c(90, 160),colours = RColorBrewer::brewer.pal(10,"RdYlGn"), labels =    mylabels, na.value = 'transparent', guide = "legend", breaks = seq(90,160,10))  
if (country == "IT"){ mapmin = mapmin +coord_fixed(xlim = c(8,10))}
mapmin = mapmin + theme_bw() + labs(title = "Estimated dates of Sowing", x = "Longitude", y = "Latitude") +
  theme(plot.title = element_text(size = 14, vjust = 1)) + labs(x = "Longitude") +
  theme(axis.text.x  = element_text(size = 8), axis.text.y  = element_text(size = 8))+
  theme(legend.position="right")+theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), 'cm'))+
  theme(legend.background = element_rect(colour = 'black'))+
  theme(legend.justification = c(0,0) ,legend.background = element_rect(colour = 'black'))
print(mapmin)

2003-2016 Sowing dates anomalies with respect to average

if (cc == 'IT'){
  gc <- geocode("pavia, it")
  center <- as.numeric(gc)
  xlims = c(8,10)
  ylims = c(44.7,45.8)
  zoom = 8
}

if (cc == 'ES'){
  gc <- geocode("castellon, es")
  center <- as.numeric(gc)
  xlims = c(-0.5,1)
  ylims = c(39,40.8)
  zoom = 8
}

if (cc == 'GR'){
  gc <- geocode("thesasloniki, gr")
  center <- as.numeric(gc)
  xlims = c(22.4, 23.55)
  ylims = c(40.4,41.2)
  zoom = 9
}

G <- ggmap(get_map(location = center,  zoom = zoom, scale = 2,maptype = 'terrain', source = 'google')) + xlim(xlims) + ylim(ylims)
mapmin = G + geom_tile(data =  newdf_avgs,aes(x = lat, y = lon, fill = anomaly), width = 0.0254 ,height = 0.0254)
mapmin <- mapmin + scale_fill_gradientn('Number of Days', limits=c(-35, 35),colours = RColorBrewer::brewer.pal(5,"RdYlGn"), na.value = 'transparent', guide = "legend", breaks = seq(-35,35,10))  
if (country == "IT"){ mapmin = mapmin +coord_fixed(xlim = c(8,10))}
mapmin = mapmin + facet_wrap(~year)
mapmin = mapmin + theme_bw() + labs(title = "Anomaly in Date of Sowing", x = "Longitude", y = "Latitude") +
  theme(plot.title = element_text(size = 14, vjust = 1)) + labs(x = "Longitude") +
  theme(axis.text.x  = element_text(size = 8), axis.text.y  = element_text(size = 8))+
  theme(legend.position="right")+theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), 'cm'))+
  theme(legend.background = element_rect(colour = 'black'))+
  theme(legend.justification = c(0,0) ,legend.background = element_rect(colour = 'black'))
print(mapmin)

2016 Sowing dates anomalies with respect to average

if (cc == 'IT'){
  gc <- geocode("pavia, it")
  center <- as.numeric(gc)
  xlims = c(8,10)
  ylims = c(44.7,45.8)
  zoom = 8
}

if (cc == 'ES'){
  gc <- geocode("castellon, es")
  center <- as.numeric(gc)
  xlims = c(-0.5,1)
  ylims = c(39,40.8)
  zoom = 8
}

if (cc == 'GR'){
  gc <- geocode("thesasloniki, gr")
  center <- as.numeric(gc)
  xlims = c(22.4, 23.55)
  ylims = c(40.4,41.2)
  zoom = 9
}

newdf_avgsub =  droplevels(subset(newdf_avgs,year %in% 2016 ))

center <- as.numeric(gc)
G <- ggmap(get_map(location = center,  zoom = zoom, scale = 2,maptype = 'terrain', source = 'google')) + xlim(xlims) + ylim(ylims)
mapmin = G + geom_tile(data =  newdf_avgsub,aes(x = lat, y = lon, fill = anomaly), width = 0.0254 ,height = 0.0254)
mapmin <- mapmin + scale_fill_gradientn('Number of Days', limits=c(-35, 35),colours = RColorBrewer::brewer.pal(5,"RdYlGn"), na.value = 'transparent', guide = "legend", breaks = seq(-35,35,10))  
if (country == "IT"){ mapmin = mapmin +coord_fixed(xlim = c(8,10))}
mapmin = mapmin + facet_wrap(~year)
mapmin = mapmin + theme_bw() + labs(title = "Anomaly in Date of Sowing", x = "Longitude", y = "Latitude") +
  theme(plot.title = element_text(size = 14, vjust = 1)) + labs(x = "Longitude") +
  theme(axis.text.x  = element_text(size = 8), axis.text.y  = element_text(size = 8))+
  theme(legend.position="right")+theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), 'cm'))+
  theme(legend.background = element_rect(colour = 'black'))+
  theme(legend.justification = c(0,0) ,legend.background = element_rect(colour = 'black'))
print(mapmin)

Data aggregated on Municipalities

The following maps show results similar to the previous ones, but aggregated on administrative units.

2003-2016 Sowing dates maps, on municipalities

if (cc == 'IT'){
  gc <- geocode("pavia, it")
  center <- as.numeric(gc)
  xlims = c(8,10)
  ylims = c(44.7,45.8)
  zoom = 8
}

if (cc == 'ES'){
  gc <- geocode("castellon, es")
  center <- as.numeric(gc)
  xlims = c(-0.5,1)
  ylims = c(39,40.8)
  zoom = 8
}

if (cc == 'GR'){
  gc <- geocode("thesasloniki, gr")
  center <- as.numeric(gc)
  xlims = c(22.4, 23.55)
  ylims = c(40.4,41.2)
  zoom = 9
}
shapefile = str_replace_all("/home/lb/projects/ermes/datasets/rs_products/Phenology/%cc%/Outputs/ERMES_Grid/Statistics/Shapefiles/ERMES_Rice_Statistics.shp","%cc%",country)
  print(shapefile)
  data_shape= readshape(shapefile)
  data_shape = spTransform(data_shape, geo_WGS84_crs)
  data_shape_df = fortify(data_shape, region = "adm_id")
  names(data_shape_df)[6] = 'adm_id'
  shape.df = join(data_shape_df, data_shape@data, by="adm_id")
  if (country != "GR") {
    shape.df_sub = droplevels(subset(shape.df, RiceFC > 10 & year %in% selected))
  } else {
    shape.df_sub = droplevels(subset(shape.df, year %in% selected))
  }
  G <- ggmap(get_map(location = center,  zoom = zoom, scale = 2,maptype = 'terrain', source = 'google')) + xlim(xlims) + ylim(ylims)
  p =  G + geom_polygon(data = shape.df_sub, aes(long,lat,group=group,fill=avgsow),color = 'grey50')
  p  <- p  + scale_fill_gradientn('Date of Sowing', limits=c(90, 160),colours = RColorBrewer::brewer.pal(10,"RdYlGn"), 
      labels = mylabels, na.value = 'transparent', guide = "legend", breaks = seq(90,160,10))  
if (country == "IT"){ p = p +coord_fixed(xlim = c(8,10))}
  p = p + facet_wrap(~year)
  p = p + theme_bw() + labs(title = "Date of Sowing", x = "Longitude", y = "Latitude") +
  theme(plot.title = element_text(size = 14, vjust = 1)) + labs(x = "Longitude") +
  theme(axis.text.x  = element_text(size = 8), axis.text.y  = element_text(size = 8))+
  theme(legend.position="right")+theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), 'cm'))+
  theme(legend.background = element_rect(colour = 'black'))+
  theme(legend.justification = c(0,0) ,legend.background = element_rect(colour = 'black'))
  print(p)

2016 Sowing dates map, on municipalities

if (cc == 'IT'){
  gc <- geocode("pavia, it")
  center <- as.numeric(gc)
  xlims = c(8,10)
  ylims = c(44.7,45.8)
  zoom = 8
}

if (cc == 'ES'){
  gc <- geocode("castellon, es")
  center <- as.numeric(gc)
  xlims = c(-0.5,1)
  ylims = c(39,40.8)
  zoom = 8
}

if (cc == 'GR'){
  gc <- geocode("thesasloniki, gr")
  center <- as.numeric(gc)
  xlims = c(22.4, 23.55)
  ylims = c(40.4,41.2)
  zoom = 9
}


if (country != "GR") {
    prova = droplevels(subset(shape.df, RiceFC > 10 & year == 2016))
  } else {
    prova = droplevels(subset(shape.df, year == 2016))
  }

  G <- ggmap(get_map(location = center,  zoom = zoom, scale = 2,maptype = 'terrain', source = 'google')) + xlim(xlims) + ylim(ylims)
  p =  G + geom_polygon(data = prova, aes(long,lat,group=group,fill=avgsow),color = 'grey50')
  p  <- p  + scale_fill_gradientn('Date of Sowing', limits=c(90, 160),colours = RColorBrewer::brewer.pal(10,"RdYlGn"), 
      labels = mylabels, na.value = 'transparent', guide = "legend", breaks = seq(90,160,10))  
if (country == "IT"){ p = p +coord_fixed(xlim = c(8,10))}

p = p + theme_bw() + labs(title = "Date of Sowing", x = "Longitude", y = "Latitude") +
  theme(plot.title = element_text(size = 14, vjust = 1)) + labs(x = "Longitude") +
  theme(axis.text.x  = element_text(size = 8), axis.text.y  = element_text(size = 8))+
  theme(legend.position="right")+theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), 'cm'))+
  theme(legend.background = element_rect(colour = 'black'))+
  theme(legend.justification = c(0,0) ,legend.background = element_rect(colour = 'black'))
print(p)

2003-2016 Sowing dates anomalies with respect to average, on municipalities

  # data_shape= readshape("/home/lb/projects/ermes/datasets/rs_products/Phenology/IT/Outputs/ERMES_Grid/Statistics/Shapefiles/ERMES_Rice_Statistics.shp")
  # data_shape = spTransform(data_in, geo_WGS84_crs)
  # data_shape_df = fortify(data_shape, region = "adm_id")
  # names(data_shape_df)[6] = 'adm_id'
  # shape.df = join(data_shape_df, data_shape@data, by="adm_id")
  # shape.df_sub = droplevels(subset(shape.df, RiceFC > 10 & year %in% selected))
  G <- ggmap(get_map(location = center,  zoom = zoom, scale = 2,maptype = 'terrain', source = 'google')) + xlim(xlims) + ylim(ylims)
  p =  G + geom_polygon(data = shape.df_sub, aes(long,lat,group=group,fill=anomsow),color = 'grey50')
  p  <- p  + scale_fill_gradientn('Number of Days', limits=c(-35, 35),colours = RColorBrewer::brewer.pal(5,"RdYlGn"), na.value = 'transparent', guide = "legend", breaks = seq(-35,35,10))  
if (country == "IT"){ p = p +coord_fixed(xlim = c(8,10))}
  p = p + facet_wrap(~year)
  p = p + theme_bw() + labs(title = "Date of Sowing", x = "Longitude", y = "Latitude") +
  theme(plot.title = element_text(size = 14, vjust = 1)) + labs(x = "Longitude") +
  theme(axis.text.x  = element_text(size = 8), axis.text.y  = element_text(size = 8))+
  theme(legend.position="right")+theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), 'cm'))+
  theme(legend.background = element_rect(colour = 'black'))+
  theme(legend.justification = c(0,0) ,legend.background = element_rect(colour = 'black'))
print(p)

2016 Sowing dates anomalies with respect to average, on municipalities

  # prova = droplevels(subset(shape.df_sub, year == 2016 & RiceFC > 10))
  G <- ggmap(get_map(location = center,  zoom = zoom, scale = 2,maptype = 'terrain', source = 'google')) + xlim(xlims) + ylim(ylims)
  p =  G + geom_polygon(data = prova, aes(long,lat,group=group,fill=anomsow),color = 'grey50')
  p  <- p  + scale_fill_gradientn('Number of Days', limits=c(-35, 35),colours = RColorBrewer::brewer.pal(5,"RdYlGn"), na.value = 'transparent', guide = "legend", breaks = seq(-35,35,10))  
if (country == "IT"){ p = p +coord_fixed(xlim = c(8,10))}

p = p + theme_bw() + labs(title = "Date of Sowing", x = "Longitude", y = "Latitude") +
  theme(plot.title = element_text(size = 14, vjust = 1)) + labs(x = "Longitude") +
  theme(axis.text.x  = element_text(size = 8), axis.text.y  = element_text(size = 8))+
  theme(legend.position="right")+theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), 'cm'))+
  theme(legend.background = element_rect(colour = 'black'))+
  theme(legend.justification = c(0,0) ,legend.background = element_rect(colour = 'black'))
print(p)