scenario Fournaise {
	
  println("Model Fournaise ready to run")
  
// Init

  // Instantiate a relation
  let rplotAg = new IsPartOf_Ag
  
  // Instantiate Study area
  let studyarea = new StudyArea
  
  // read Plots from Shapefile  
  let plotShpAg = new ShpAg
  let lplotAg = plotShpAg.readAll()  // returns a list of PlotAg
  for (p:lplotAg) {
    	rplotAg.connect(p, studyarea) 
  }
  
  // Instantiate neighbourhood graph and connect those plots less than 20m apart
  let gN = new Neighbourhood
  gN.addAllPlotAg(lplotAg)
  gN.complete.near(20).connect()
  

  // Initialize kml file for output
  let kmln = new Kml_LC
  let plt = colorRange(8,"dem2")        // get list of colors from predefined palette "dem2"
  kmln.defStyleRange("st",1.0,plt,-0.2) // prefix, line thickness, list of colors , darken line color a little
  kmln.defStyle("neighb",1.5,Color|rgb(66,149,217),Color|rgb(66,149,217)) // set color for neighbourhood graph
  
  // Initialize shapefile for output
  let shpout = new ShpAgOut
  
// End Init



// Simulation
   // 
   let ybeg = 2014
   let yend = 2030
   for (year : ybeg..yend) {
      println(year)
      
      //Initialise list of plots to be added to output shapefile 
      let lplotAgOut = new List<PlotAg>

      // Write land cover state in kml file and shapefile
      for (p:lplotAg) {
    	  p.outputKml(year, kmln) // write to Kml for each PlotAg of the list
    	  
    	  p.date_begin = ""+ year +"-01-01"
    	  p.date_end = ""+ year +"-12-31"
    	  lplotAgOut.add(p)      // add plot to list of plots for shapefile output
   	  
      }
      // output to shapefile
      shpout.append(lplotAgOut)

      
      // Change land cover according to neighbourhood
      gN.changeLC()

             
   }  


// Saving kml file 
  println("")
  print("Saving kml file...")
  gN.drawEdge(ybeg,yend,kmln)
  kmln.saveAsKml()
  println("ok.")	

// Saving shapefile  
  print("Saving shp ... ")
  shpout.close()
  println("ok.")

  println("Done.")
  
  
}



// Define the agricultural plot entity. 
entity PlotAg {
	property MultiPolygon geom
	property Integer id
	property Integer lct // land cover type
	property Integer surf // surface area in mē
	property String date_begin // for use with time manager plugin
	property String date_end  //  in QGIS

	service outputKml (Integer year, KmlExport kml){
        let id = "Plot" + id + "_" + year
        let style="st" + lct 
        let deb = year+"-01-01"
        let fin = ""+(year+1)+"-01-01"//"-12-31"
        kml.addGeometry (style,id, deb, fin, geom, style,0) 
        
    }
    

	
}

// Define the study area entity. 
entity StudyArea {
	
 }


// Define a datafacer to read shapefile of plots 
datafacer ShpAg {
	data Shapefile("data/AgPlots_Fournaise.shp","EPSG:32740") // WGS84 UTM40S

	// 
	match PlotAg { 
  	              geom : "geom"
  	              id : "PID"
   	              lct : "CCULTURE"
   	              surf : "SURF"
   	                     	          				
	}
}

datafacer ShpAgOut {
	data Shapefile("output/Fournaise_Landcover.shp","EPSG:32740")
	match PlotAg {
		          geom : "the_geom"
		          id : "gid"
		          lct : "lct"
		          date_begin : "begin"
		          date_end : "end"
	}
}



// Define a datafacer for Kml export kml - Land cover
datafacer Kml_LC {
	data KmlExport("output/Fournaise_LandCover.kml")
}



// Define relation between agricultural plots and study area

relation IsPartOf_Ag <PlotAg p, StudyArea s> {
	
	interaction outputKml (Integer year, KmlExport kml){
        let id = "Plot" + p.id + "_" + year
        let style="st" + p.lct 
        let beg = year+"-01-01"
        let end = ""+(year+1)+"-01-01"//"-12-31"
        kml.addGeometry (style,id, beg, end, p.geom, style,0) 
        
    }	
}
	
// aggregation operator for changing the lct of a plot according to that of its neighbours
aggreg Integer Majority {	
	let nbn = values.size()
	let f0 = 0
	let maj = 0
	if (nbn > 0) {
		 // select the most frequent lct among neighbours
	     for (i:1..8) {
	           let f = values.frequency(i)
	           if (f > f0) {
	   	           f0 = f
	   	           maj = i
	           }
	     }	
	} 
	// but if more than half of the neighbours have the same lct, it takes a random one between 1 and 8
	let r = 1.0 + random()*7.5
	if ((f0 > nbn*0.5) && (random() > 0.2)){
		maj = r.intValue
	}
		
	return maj
    }
	
relation Neighbourhood<PlotAg p1, PlotAg p2>{

	filter near(Integer dist){
		return (p1.geom.envelope.distance(p2.geom.envelope) < dist) &&
		 (p1.geom.distance(p2.geom) < dist)
	}
	
	interaction changeLC() {
		// the land cover type of a plot is changed according to the lct of its neighbours
		p1.lct =  p2.lct  // a plot receives the lct of the neighbouring plot
		p2.lct =  p1.lct 
		
	} agg {
		p1.lct << Majority  // a plot has many neighbours, therefore an aggregation function is applied 
		p2.lct << Majority  // default functions (Mean, Min, Max, Sum ...) or user defined (here Majority)
	}
	
	
	interaction drawEdge(Integer yearbeg, Integer yearend, KmlExport kml){
        let pc1 = Point|xy(p1.geom.centroid.x,p1.geom.centroid.y)
        let pc2 = Point|xy(p2.geom.centroid.x,p2.geom.centroid.y)
		let arc = Line|points(pc1,pc2)
		let beg = yearbeg+"-01-01"
        let end = ""+(yearend+1)+"-01-01"//"-12-31"
		kml.addGeometry("Neighbours",p1.id+" -- "+p2.id,beg,end,arc,"neighb",0)
    }
}