rm(list = ls(all = TRUE))
library(akima)
library(maps)
library(mapproj)
library(fields)
library(mapdata)
library(maptools)
library(spam)
library(raster)
library(plotrix)
library(mapplots)
library(ggmap)
source('INPUT_R/function.R')
library(dismo)
library(colorRamps)


#fichier avec les frontieres
#shp<-readShapeSpatial("INPUT_R/ntlbndscarto")

#exemple lon lat data
data<-read.table("exemple",sep=',',header=T)
#station,lon,lat,R
#A1,2.35652,48.92476,0.519937912172249
#AUT,2.2533,48.8503,0.552200410830095
#BASCH,2.3275,48.8278,0.530422476320368
lon<-data$lon
lat<-data$lat
var<-data[,4]


#definir le cadre
lonmin<-min(lon)
lonmax<-max(lon)
latmin<-min(lat)
latmax<-max(lat)

#recupération du fond de carte
e <- extent(lonmin-1 , lonmax+1 , latmin-0.1 , latmax+0.1)  ## extenhelpt (long, long, lat, lat)
r <- gmap(e, type="satellite", scale=2,lonlat = TRUE)


nline<-length(lon)
#plot
heights=c(5)
widths=c(5,1.2)
res=200
zlim<-c(0,round(max(data[,4],na.rm=T)))
legendlab<-colnames(data)[4]
titlelab<-colnames(data)[4]
png(filename = "map_points.png", res=res,width = sum(widths), height = sum(heights), units="in",pointsize=18)
   layout(matrix(c(1,2),nrow=1,ncol=2,byrow=TRUE), widths = widths, heights = heights, respect=TRUE)
   plot(c(lonmin,lonmax),c(latmin,latmax),type='n',xlab="Longitude",ylab="Latitude")
   #plot google map
   plot(r,add=T)
   #plot boundaries
   #plot(shp,add=T)
   map(add=T)
   map('france',add=T)
   #convert val to color
   my_col<-colorRamps::matlab.like(100)
   col<-val2col(var,col=my_col,zlim=zlim)
   #plot points
   for (i in seq(1,nline)){
       if(is.na(data[i,4]))next
       points(lon[i],lat[i],col="black",pch=21,bg=col[i],cex=1)
   }
   title(titlelab)
   par(mai=c(0.4, 0, 0.4, 1))
   image.scale(var,zlim=zlim, horiz=FALSE, yaxt="n",col=my_col)
   axis(4,las=1)
   mtext(legendlab, side=4, line=2.5)
   box()
dev.off()