// DATA PREPARATION
//Load the nightlight data
var light_dataset = ee.ImageCollection('NOAA/DMSP-OLS/NIGHTTIME_LIGHTS').select('stable_lights');
Map.addLayer(light_dataset,{min: 3.0,max: 60.0,palette:['000000','9999CC']},'nighttimeLights');
//Load province boundary of China
var China = China;
Map.addLayer(China,{color:'FFFFFF'},'China_Province');
Map.centerObject(China,4);
print(China);
//Load country boundary of China
var boundary_geo=China.geometry();
var boundary = ee.Feature(boundary_geo);


// Part.1 AN OVERVIEW OF URBAN BUILT-UP AREA CHANGES IN CHINA
//Load the nightlight from 2003-2013
var light_2003 = ee.Image('NOAA/DMSP-OLS/NIGHTTIME_LIGHTS/F152003').select('stable_lights').clip(boundary);
Map.addLayer(light_2003,{min: 3.0,max: 60.0,palette:['000000','FF9933','ffffff'],opacity:0.8},'Lights_2003');
var light_2008 = ee.Image('NOAA/DMSP-OLS/NIGHTTIME_LIGHTS/F162008').select('stable_lights').clip(boundary);
Map.addLayer(light_2008,{min: 3.0,max: 60.0,palette:['000000','FF9933','ffffff'],opacity:0.8},'Lights_2008');
var light_2013 = ee.Image('NOAA/DMSP-OLS/NIGHTTIME_LIGHTS/F182013').select('stable_lights').clip(boundary);
Map.addLayer(light_2013,{min: 3.0,max: 60.0,palette:['000000','FF9933','ffffff'],opacity:0.8},'Lights_2013');

// find the place that are considered "lit"
var blank = ee.Image(0);
function AreaLit(i){
  var threhold=blank.where(i.gt(31),1);
  var alit = threhold.mask(threhold);
  alit = ee.Image.pixelArea().mask(alit);
  alit=alit.divide(1000000);
  alit=alit.set('index',i.get('system:index'));
  return alit;
}
// stable nighttime light area change (by province) from 2003-2013
function tabulate(i){
  return China.map(function(f){
    var r=i.reduceRegion({
      reducer:ee.Reducer.sum(),
      geometry:f.geometry(),
      scale:500,
      bestEffort:true,
      maxPixels:1e9});
    return ee.Feature(null,{
      name:f.get('Province'),
      area:r.get('area'),
      index:i.get('index')});});}
var arelit=light_dataset.map(AreaLit);
var arelit_table=arelit.map(tabulate).flatten();
Export.table.toDrive(arelit_table,'DMSP_results','DMSP results','results','csv');


// Part.2 CALCULATE THE POSITIVE/NEGATIVE GROWTH OF URBAN BUILT-UP AREAS
// Add a band containing image data as years
function createTimeBand(img){
  return img.addBands(img.metadata('system:time_start').divide(1e18));}
//fit a linear trend to the nighttime lights collection
var getyear=function(year){
  return ee.Date.fromYMD({
    day:1,month:1,year:year});};
var linear_fit=function(y1,y2){
  var startyear = getyear(y1);
  var endyear = getyear(y2);
  var timeband=light_dataset.filterDate(startyear,endyear).map(createTimeBand);
  var linearfit=timeband.select(['system:time_start','stable_lights']).reduce(ee.Reducer.linearFit()).clip(boundary);
  return linearfit;
};
var fit_china_92_95=linear_fit(1992,1995);
var fit_china_96_10=linear_fit(1996,2010);
var fit_china_10_14=linear_fit(2010,2014);
// Display Trend in Red(positive growth),Blue(negative growth), and green(stable brightness)
var visParams={min:0,max:[0.18,20,-0.18],bands:['scale','offset','scale'],opacity:0.8};
Map.addLayer(fit_china_92_95,visParams,"linear_fit_92_95");
Map.addLayer(fit_china_96_10,visParams,"linear_fit_96_10");
Map.addLayer(fit_china_10_14,visParams,"linear_fit_10_14");

// Visualize the rate of urban sprawl and shrinkage by province
// extract the pixels with positive growth (The slope of linear regression is positive) and negative growth (The slope of linear regression is negative)
var positive_growth=function(img){
  var img_scale=img.select("scale");
  var positive=img_scale.multiply(img_scale.gt(0));
  return positive;
};
var negative_growth=function(img){
  var img_scale=img.select("scale");
  var negative=img_scale.multiply(img_scale.lt(0).multiply(-1));
  return negative;
};
var positive_10_14=positive_growth(fit_china_10_14);
var negative_10_14=negative_growth(fit_china_10_14);
Map.addLayer(positive_10_14,{min: 3.0,max: 60.0,palette:['000000','FF6666'],opacity:0.8},"positive_10_14");
Map.addLayer(negative_10_14,{min: 3.0,max: 60.0,palette:['000000','006699'],opacity:0.8},"negative_10_14");

//calculate the area of positive_growth and negative_growth by province, respectively
function count_positive(Chi){
  var area =positive_10_14.reduceRegion({
    reducer:ee.Reducer.sum(),
    geometry:Chi.geometry(),
    scale:500,
    bestEffort:true,
    maxPixels:1e9});
  var Chi_copy=Chi;
  return Chi_copy.set({LitArea:area.get('scale')});
}
function count_negative(Chi){
  var area =negative_10_14.reduceRegion({
    reducer:ee.Reducer.sum(),
    geometry:Chi.geometry(),
    scale:500,
    bestEffort:true,
    maxPixels:1e9});
  var Chi_copy=Chi;
  return Chi_copy.set({LitArea:area.get('scale')});
}
var area_positive=China.map(count_positive);
var area_negative=China.map(count_negative);
//calculate the ratio of positive_growth/negative_growth area and total area
function density(Chi){
  var Area=ee.Number(Chi.area());
  var ratio=ee.Number(Chi.get('LitArea')).divide(Area);
  var Chi_copy=Chi;
  return Chi_copy.set({"ratio":ratio});}
var ratio_positive=area_positive.map(density);
var ratio_negative=area_negative.map(density);
// plot the choropleth map of ratio by province
var choropleth=function(fc){
  var img=fc.reduceToImage(["ratio"],ee.Reducer.first());
  return img;
};
var img_positive=choropleth(ratio_positive);
Map.addLayer(img_positive,{max:130,min:0,palette:["FFFFFF","E99C9A","E96D6A","C02A26","5D0606"],opacity:1},"positive growth area/total area");
var img_negative=choropleth(ratio_negative);
Map.addLayer(img_negative,{max:130,min:0,palette:["FFFFFF","8FB2D1","679FD1","20517C","092C4C"],opacity:1},"negative growth area/total area");
//identify provinces with the largest urban expansion rate and contraction rate
var top_positive=ratio_positive.limit(10,"ratio",false);
print(top_positive);
var top_negative=ratio_negative.limit(10,"ratio",false);
print(top_negative);
//export the result table
Export.table.toDrive({
  collection:top_positive,
  description: "top_positive",
  fileNamePrefix: "top_positive",
  fileFormat: "KML",
  selectors: ["Province","ratio"]
});
Export.table.toDrive({
  collection:top_negative,
  description: "top_negative",
  fileNamePrefix: "top_negative",
  fileFormat: "KML",
  selectors: ["Province","ratio"]
});

// Part.3 CHARACTERISTIC OF URBAN CONTRACTION & EXPANSION AREAS
//Load the selected two province
var selected = selected;
Map.addLayer(selected,{color:'FFFFFF'},'selected_province');
Map.centerObject(selected,4);
print(selected);