<jupyter_start><jupyter_text>**Irrigated Area Classification in Semi-Arid Regions**This notebook contains an example of how to classify irrigated herbaceous and perennial tree areas in semi-arid regions. It makes use of the Google Earth Engine, for which an account is needed. For more information on the Google Earth Engine and instructions on how to make an account click [here](https://earthengine.google.com/).   First, Import the Github Repository containing the Irrigated Areas Classification package.<jupyter_code>!git clone https://github.com/Remdeht/ia_detector.git ia
 !pip install monthdelta<jupyter_output><empty_output><jupyter_text>Load the required libraries and initialize the GEE. If this is your first tume using the Earth Engine, you need to authenticate first via the *ee.Authenticate()* call. Instructions on how to authenticate your GEE will appear in the output below.<jupyter_code>import ee 
ee.Authenticate()  # if first time user of the Google Earth Engine on this device, otherwise ee.Initialize() is enough
import ee
ee.Initialize()

import itertools
import sys
import numpy as np
import folium
from ia.gee_functions.classification import create_features, create_training_areas, classify_irrigated_areas, \
    join_seasonal_irrigated_areas
from ia.gee_functions.constants import GEE_USER_PATH
from ia.gee_functions.export import track_task, export_to_drive
from ia.gee_functions import visualization
from datetime import datetime
from ia.gee_functions.validation import calc_area<jupyter_output><empty_output><jupyter_text>**Area of Interest**Before running the classification you must select an area of interest. You can use your own area, by filling in the Lat/Long coordinates of the rectangular extent of your choice or use one of the extents available from the list below. If you run the next code block a map will appear that allows you to check the coordinates of a location by clicking on it.<jupyter_code>map = visualization.create_folium_map(images=None, coords=[37.656830, -1.046193], zoom=8, height='100%')
pu = folium.LatLngPopup().add_to(map)
map<jupyter_output><empty_output><jupyter_text>In the dropdown menu, select which area of interest you'd wish to use for the classification. If you'd like define your own extent, select the option 'Custom Rectangular Extent'. You can then insert the maximum and minimum latitude and longitude coordinates which will be used to create a rectangular polygon.If you're familiar with the Google Earth Engine, and know how to upload vector layers and load them via their assetId, feel free to use your own vector asset by selecting the 'Other GEE Vector' option in the drop down menu and specifying the assetId in the 'aoi_loc' field.<jupyter_code>#@markdown ---
#@title Select your area of interest { run: "auto" }
area = "Campo de Cartagena" #@param ["Campo de Cartagena", "Murcia Region", "Segura Watershed", "Custom Rectangular Extent", "Other GEE Vector"] { run: "auto" }

areas = {
    'Campo de Cartagena':'cdc',
    'Murcia Region':'rdm',
    'Segura Watershed':'cds',
    'Custom Rectangular Extent':'',
    'Other GEE Vector':''
}

area_of_interest = areas[area]

if not area_of_interest == '':
  aoi = ee.FeatureCollection(f'users/Postm087/vector/outline/outline_{area_of_interest}')

#@markdown ---
#@markdown ### *Custom Rectangular Extent*
#@markdown ### Fill in maximum and minimum latitude and longitude coordinates of the extent
elif area == 'Custom Rectangular Extent':
  lat_min = 37.656830  #@param {type:"number"}
  lat_max = 37.784704  #@param {type:"number"}
  lng_min = -1.046193 #@param {type:"number"}
  lng_max = -0.832859 #@param {type:"number"}
  aoi = ee.FeatureCollection(ee.Feature(ee.Geometry.Polygon([
                   [lng_min, lat_max],
                   [lng_max, lat_max],
                   [lng_max, lat_min],
                   [lng_min, lat_min],
                   [lng_min, lat_max]
  ])))
  #@markdown ### Fill in the name of your area, this will be used to name the classification outputs { run: "auto" }
  area_of_interest = 'cre' #@param {type:"string"}
#@markdown ---
#@markdown ### *Load in your own Vector via its assetId*
elif area == 'Other GEE Vector':
  aoi_loc = ''#@param {type:"string"}
  aoi = ee.FeatureCollection(aoi_loc)
  area_of_interest = 'cre' #@param {type:"string"}


aoi_coordinates = aoi.geometry().bounds().getInfo()['coordinates']
aoi_centroid = aoi.geometry().bounds().centroid(1).coordinates().getInfo()

layer = {
    'Area of Interest':aoi.getMapId({'color':'White'})
}

map = visualization.create_folium_map(layer, coords=[aoi_centroid[1], aoi_centroid[0]], zoom=8, height='100%')
map<jupyter_output><empty_output><jupyter_text>Now that the area for classification has been specified, the next step before classification is to determine the period of time for which to perform the classification. Two separate classifications are performed, one for the summer season which spans from April until September and once for the winter which spans from October until March. You can specify the period of your choosing in the date-picker below, but keep in mind to select a one year period.<jupyter_code>start_date = "2013-04-01" #@param {type:"date"}
end_date = "2014-04-01" #@param {type:"date"}


start_date_dt = datetime.strptime(start_date, '%Y-%m-%d')
end_date_dt = datetime.strptime(end_date, '%Y-%m-%d')
year = str(end_date_dt.year)[-4:]
classification_period = {year:(start_date, end_date)}

# Some variables needed for classification but not of much importance for the user.
sat = 'landsat'  # Satellite data to use, as of now sentinel has not been tested properly so better not to use it. 
stats = ['mean', 'max'] # Statistical maps to use for classification.
stats_combos = list(itertools.combinations(stats, 2))

# Some folders where the data are saved
crop_data_folder = f'{GEE_USER_PATH}/ia_classification/raster/data/{area_of_interest}/{sat}/'
training_data_folder = f'{GEE_USER_PATH}/ia_classification/raster/training_areas/{area_of_interest}/'
results_folder = f'{GEE_USER_PATH}/ia_classification/raster/results/'<jupyter_output><empty_output><jupyter_text>Now that both the area and the time period for classification are known, the feature data for classification can be generated from satellite imagery. For both the summer and winter season a collection of pixel statistic maps will be generated based on a selection of spectral indices. The spectral indices used for classification are:1.   Normalized Difference Vegetation Index - *NDVI*2.   Normalized Difference Water Content Index - *NDWI*2.   Normalized Difference Water Bodies Index - *NDWI*3.   Water-Adjusted Green Index - *WGI*4.   Green Chlorophile Vegetation Index - *GCVI*5.   Normalized Difference Built-Up Index - *NDBI*The generation of the feature data maps can take a while, depending on the size of the area of interest. When the task is done running on the Google Earth Engine Servers, a message will pop up in the output below.<jupyter_code>tasks = {}
for year in classification_period:
    export_feature_task = create_features(
        classification_period[year],
        aoi,
        aoi_name=area_of_interest,
        sensor=sat)
        
track_task(export_feature_task) # tracks the status of the export task<jupyter_output><empty_output><jupyter_text>Let's have a look at some of the feature data layers. By running the code block below you'll generate a folium map containing the RGB composites for the summer and the winter seasons, as well as the median NDVI values and the WGI standard deviation values per pixel over the summer and winter.<jupyter_code>summer = ee.Image(f"{crop_data_folder}crop_data_summer_{area_of_interest}_{year}")
winter = ee.Image(f"{crop_data_folder}crop_data_winter_{area_of_interest}_{year}")
images = {
      'Summer Feature Data RGB': summer.getMapId(visualization.vis_params_rgb(bands=['red', 'green', 'blue'])), 
      'Winter Feature Data RGB': winter.getMapId(visualization.vis_params_rgb(bands=['red', 'green', 'blue'])),
      'Summer Feature Data NDVI mean': summer.getMapId(visualization.vis_params_cp(band=['NDVI_mean'], min_val=-1, max_val=1)), 
      'Winter Feature Data NDVI mean': winter.getMapId(visualization.vis_params_cp(band=['NDVI_mean'], min_val=-1, max_val=1)), 
      'Summer Feature Data WGI standard deviation': summer.getMapId(visualization.vis_params_cp(band=['WGI_std'], min_val=0, max_val=1)), 
      'Winter Feature Data WGI standard deviation': winter.getMapId(visualization.vis_params_cp(band=['WGI_std'], min_val=0, max_val=1)), 
  }
map = visualization.create_folium_map(images, coords=[aoi_centroid[1], aoi_centroid[0]], zoom=10, height='100%')
map<jupyter_output><empty_output><jupyter_text>**Training Areas**Now that the feature data is generated, the next step in to create training areas for the Random Forest thresholding. To save you the task of manually selecting these areas, a system of automatically selecting training areas for classification has been applied. The system makes use of a priori knowledge of typical spectral characteristics for a selection of land cover classes to select training sites via thresholding.Training sites for the following land cover classes will be generated: *   Forest*   Scrubs*   Rainfed Crops and Trees*   Greenhouses*   Irrigated Herbaceous Crops*   Irrigated Perennial Trees*   Water Bodies*   Unproductive Areas (Urban/Fallow Lands)By running the code block below the generation of training sites will be started. This may take a couple of minutes. A confirmation message will appear when the process is completed.<jupyter_code>training_tasks = {}
for year in classification_period:
    export_training_areas_task = create_training_areas(
          aoi,
          f"{crop_data_folder}crop_data_season_{area_of_interest}_{year}",
          aoi_name=area_of_interest,
          year_string=year
          )
track_task(export_training_areas_task)  # tracks the export task<jupyter_output><empty_output><jupyter_text>If you've reached this point, the generation of training areas should have been completed. Run the code block below to see the results.<jupyter_code>training_classes = {
    0:'Unused',
    1:'Forest',
    2:'Scrubs',
    3:'Rainfed Crops and Trees',
    4:'Greenhouses',
    5:'Irrigated Herbaceous Crops',
    6:'Irrigated Trees',
    7:'Water Bodies',
    8:'Urban/Fallow lands'
}

summer_training = ee.Image(f'{training_data_folder}training_areas_summer_{area_of_interest}_{year}')
winter_training = ee.Image(f'{training_data_folder}training_areas_winter_{area_of_interest}_{year}')
images_training = {
      'Summer Feature Training Areas': summer_training.getMapId(visualization.vis_rf_classification(band='training')),
      'Winter Feature Training Areas': winter_training.getMapId(visualization.vis_rf_classification(band='training')) 
  }
map = visualization.create_folium_map(images_training, coords=[aoi_centroid[1], aoi_centroid[0]], zoom=10, height='100%')
visualization.create_categorical_legend(map, visualization.vis_rf_classification()['palette'], training_classes)
map<jupyter_output><empty_output><jupyter_text>**Classification**Now that the training sites and feature data are generated, we can finally move on to the classification. The classification is performed using a Random Forest classifier, for which the following parameters can be adjusted: *  The Number of Trees - *no_trees**  The Number of Variables per Split - *vps**  The Bagging Fraction - *bf**  Minimum Number of Training Points per class - *min_tp**  Maximum Number of Training Points per class - *max_tp*Run the code below to start the classification. The classification may make take some time depending on the number of trees and training points specified.<jupyter_code>#@markdown ---
#@title ### Classification Parameters
#@markdown ### Number of Trees
no_trees = 150 #@param {type:"slider", min:0, max:500, step:10}
#@markdown ### Variables per Split
vps = 6 #@param {type:"slider", min:2, max:10, step:1}
#@markdown ### Bagging Fraction
bf = 0.2 #@param {type:"slider", min:0, max:1, step:0.05}
#@markdown ### Minimum Number of Training Points
min_tp = 1000 #@param {type:"slider", min:100, max:5000, step:100}
#@markdown ### Maximum Number of Training Points
max_tp = 5000 #@param {type:"slider", min:5000, max:75000, step:1000}
classification_tasks = {}

for year in classification_period:
  for season in ['winter', 'summer']:
    for combo in stats_combos:
      # load the feature data maps
      feature_data_image = ee.Image(f"{crop_data_folder}crop_data_{season}_{area_of_interest}_{year}")
      bands_to_select = ['red', 'green', 'blue', 'nir', 'swir1','.*std.*', 'MTI']
      stat_bands = [f'.*{s}.*' for s in list(combo)]
      bands_to_select += stat_bands # bands to select for classification
      feature_data_image = feature_data_image.select(bands_to_select)
      classification_name = "_".join(combo)
      classification_name += f'_{season}'
                  
      training_image = ee.Image(
      f'{GEE_USER_PATH}/ia_classification/raster/training_areas/{area_of_interest}/training_areas_{season}_{area_of_interest}_{year}')
      classification_task = classify_irrigated_areas(
      feature_data_image,
      training_image,
      aoi,
      aoi_name=area_of_interest,
      season=season,
      year=year,
      no_trees=no_trees,
      bag_fraction=bf,
      vps=vps,
      min_tp=min_tp,
      max_tp=max_tp
      )
      classification_tasks[season]=classification_task
          
track_task(classification_tasks)<jupyter_output><empty_output><jupyter_text>Now that the classification has been completed, the results can be inspected using a folium map.<jupyter_code>summer_clf = ee.Image(f'{results_folder}random_forest/{area_of_interest}/ia_random_forest_{no_trees}tr_{vps}vps_{int(bf*100)}bf_{area_of_interest}_summer_{year}')
winter_clf = ee.Image(f'{results_folder}random_forest/{area_of_interest}/ia_random_forest_{no_trees}tr_{vps}vps_{int(bf*100)}bf_{area_of_interest}_winter_{year}')
images_clf = {
      'Summer Feature Data NDVI mean': summer.getMapId(visualization.vis_params_cp(band=['NDVI_mean'], min_val=-1, max_val=1)), 
      'Winter Feature Data NDVI mean': winter.getMapId(visualization.vis_params_cp(band=['NDVI_mean'], min_val=-1, max_val=1)),
      'Summer Feature Data WGI standard deviation': summer.getMapId(visualization.vis_params_cp(band=['WGI_std'], min_val=0, max_val=1)), 
      'Winter Feature Data WGI standard deviation': winter.getMapId(visualization.vis_params_cp(band=['WGI_std'], min_val=0, max_val=1)),
      'Summer RF Classification results': summer_clf.getMapId(visualization.vis_rf_classification()),
      'Winter RF Classification results': winter_clf.getMapId(visualization.vis_rf_classification())
}

map = visualization.create_folium_map(images_clf, coords=[aoi_centroid[1], aoi_centroid[0]], zoom=10, height='100%')
map = visualization.create_categorical_legend(map, visualization.vis_rf_classification()['palette'], training_classes)
map<jupyter_output><empty_output><jupyter_text>Finally, the last step is to combine the classification results for the summer and the winter season into a single map by running the code block below.<jupyter_code>ia_summer = summer_clf.select('irrigated_area')
ia_winter = winter_clf.select('irrigated_area')
task = join_seasonal_irrigated_areas(
    ia_summer,
    ia_winter,
    area_of_interest,
    year,
    aoi,
    export_method='asset',
  )
track_task(task)<jupyter_output><empty_output><jupyter_text>**Results**The irrigated areas from both classification maps have been combined, creating a final map depicting the irrigated areas over the whole year. The following classes are assigned, based on which period the area was classified as irrigated. *  Not Irrigated*  Year-Round Irrigated Trees*  Year-Round Irrigated Crops*  Summer Irrigated Trees*  Summer Irrigated Crops*  Winter Irrigated Trees*  Winter Irrigated Crops*  Uncertain AreasLet's run the code below and take a look at the final result!<jupyter_code>ia_classes = {
        0: 'Not Irrigated',
        1: 'Year Round Irrigated Trees',
        2: 'Year Round Irrigated Crops',
        3: 'Summer Irrigated Trees',
        4: 'Summer Irrigated Crops',
        5: 'Winter Irrigated Trees',
        6: 'Winter Irrigated Crops',
        7: 'Uncertain Areas',
    }


ia_year = ee.Image(f'{results_folder}irrigated_area/{area_of_interest}/irrigated_areas_{area_of_interest}_{year}').clip(aoi)

images_results = {
      'Irrigated Areas Overview': ia_year.getMapId(visualization.vis_irrigated_area_map()),
  }
map = visualization.create_folium_map(images_results, coords=[aoi_centroid[1], aoi_centroid[0]], zoom=10, height='100%')
map = visualization.create_categorical_legend(map, visualization.vis_irrigated_area_map()['palette'], ia_classes)
map
#@title Export Results to Drive {run:"auto"}
#@markdown If you'd like to export the classification results to your Google drive account, please check the box below and run the code. 

export_results_to_drive = True #@param {type:"boolean"}


images_to_export = {
                    'summer_classification': summer_clf,
                    'winter_classification': winter_clf,
                    'irrigated_areas': ia_year,
                    }
                  
if export_results_to_drive:
  for img in images_to_export:
    export_task = export_to_drive(images_to_export[img], 'image', img, aoi_coordinates, 'ia_classification' )
  track_task(export_task)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09821734658896127, 'fraction_numerical': 0.011598674437207177, 'mean_word_length': 5.150035137034434, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '292590', 'license': 'MIT', 'n_tokens_mistral': 5596, 'n_tokens_neox': 5266, 'n_words': 1615, 'path': 'IA_Classification_Script.ipynb', 'repo_name': 'Remdeht/ia_detector', 'size': '17502'}