<jupyter_start><jupyter_text>0. FUNCTIONS<jupyter_code># Function to add year attribute to image
def addyear(img):
    return img.set('year', ee.Image(img).date().get('year'))

# Function to add Day Of Year attribute to image
def addDOY(img):
    return img.set('DOY', ee.Image(img).date().getRelative('day', 'year'))

# Function to full date
def addDATE(img):
    return ee.Image.constant(img.date().getRelative('day', 'year')).int().updateMask(img.select(0).mask())


# Parameters for masking function (move to preamble)
def landsatmasking(img):

  qa = img.select('pixel_qa')
  # If the cloud bit (5) is set and the cloud confidence (7) is high
  # or the cloud shadow bit is set (3), then it's a bad pixel.
  cloudorsnow = qa.bitwiseAnd(1 << 5). \
    And(qa.bitwiseAnd(1 << 7)). \
    Or(qa.bitwiseAnd(1 << 3)). \
    Or(qa.bitwiseAnd(1 << 4))  
  # Remove edge pixels that don't occur in all bands
  mask2 = img.mask().reduce(ee.Reducer.min())
  return img.updateMask(cloudorsnow.Not()).updateMask(mask2)

# Spectral indices function: adds index as a new band to every image
def addINDICES(img):
    
    ndvi = img.normalizedDifference(['B4', 'B3']).rename('NDVI')
    lswi = img.normalizedDifference(['B4', 'B5']).rename('LSWI')
    ndsi = img.normalizedDifference(['B2', 'B5']).rename('NDSI')    
    # I multiply by scale factor for EVI (DIDNT WORK TO SOLVE HIGH VALUES ISSUE)
    evi = img.expression(
    '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
        'NIR': img.select('B4').multiply(0.0001),
      'RED': img.select('B3').multiply(0.0001),
      'BLUE': img.select('B1').multiply(0.0001)
    }).rename('EVI')
    
    return img.addBands(ndvi).addBands(lswi).addBands(ndsi).addBands(evi)

# Function to add a band of rice classification to each image in the col
# Rule for rice is that LSWI − NDVI > 0 or LSWI − EVI > 0
def addRICE(img):    
    rice = img.expression(
    "(lswi - ndvi > 0) || (lswi - evi > 0) ? 1\
    :2",{'lswi': img.select('LSWI'),'ndvi': img.select('NDVI'),'evi': img.select('EVI')}).rename('RICE')
    
    return img.addBands(rice)

# Apply Masks Function
def landcovermasking(img):
  return img.updateMask(mask1.Not()).updateMask(mask2.Not()).updateMask(mask3.unmask().Not()).updateMask(mask4.Not())

# Function to create binary layer for rice class
def classornot(img):
    return img.eq(1)

# Percentile calculator
def percentileextract(datedif):
    if datedif < 43 :
        percentile = 70
    elif datedif > 70 :
        percentile = 90
    else :
        percentile = 80
    return percentile<jupyter_output><empty_output><jupyter_text>1. INPUTS AND TABLE LOADING<jupyter_code># Out Directory
maps_outdir = os.path.join(os.path.expanduser('~'), '../rice-mapping-japan/results/maps')

# Map type: binary or representing rice types (normal, early, second transplants'
# 'yes' for binary, 'no' for maps with different types
# If not binary, map values are: 1: normal rice; 10: early; 20: second; other values are combinations of those (rice detected in different periods)
binarymap = 'yes'  

# Define periods names
periodnames = {'8589','9094','9599','0004','0509','1014','1519'}<jupyter_output><empty_output><jupyter_text>Insert Japan Map<jupyter_code>Map = geemap.Map(center=[38,138], zoom=6)
Map.add_basemap('HYBRID')
Map<jupyter_output><empty_output><jupyter_text>Load prefectures shapefile and visualize<jupyter_code># Japan shp is not small, so we ingest it from personal account (otherwise we can upload locally with geemap.shp_to_ee)
japan_shp = 'users/luiscartor/japan_gadm0'
japan = geemap.ee.FeatureCollection(japan_shp)

point = ee.Geometry.Point([38,138])
Map.addLayer(point, {}, 'Point')

# Load prefecture shapefile (from GEE account)
prefs_shp = 'users/luiscartor/japan_gadm1_noOKINAWA'
prefs = geemap.ee.FeatureCollection(prefs_shp)

Map.addLayer(prefs, {}, 'Prefectures')

# Can use this block as playground visualize auxiliary data layers (e.g. masks) created later<jupyter_output><empty_output><jupyter_text>Data Tables  Transplanting data<jupyter_code># List of prefectures transplanting START DOY
# These data can be found in Carrasco et al 2022 (In prep.)
transplant_startdates_table = '../rice-mapping-japan/data/transplanting_startdates.csv'
transplant_startdates = pd.read_csv(transplant_startdates_table, sep=",")
transplant_startdates.head()
# Due to data inconsistencies between years, we use starting dates from 9094 for 8589
transplant_startdates['8589'] = transplant_startdates['9094']
transplant_startdates.head()
# List of prefectures transplanting START DOY
# This data can be found in Carrasco et al 2022 (In prep.)
transplant_enddates_table = '../rice-mapping-japan/data/transplanting_enddates.csv'
transplant_enddates = pd.read_csv(transplant_enddates_table, sep=",")
transplant_enddates.head()
# Trasplanting period lenghts
transplantdates_differences = pd.DataFrame({'Prefecture':transplant_enddates['Prefecture'],'difference':transplant_enddates['8589']+30 - transplant_startdates['8589']})
transplantdates_differences.head()
# Some prefectures have "early" or "second" transplantations
# Here we load the dates of those especial transplantation periods
# These data can be found in Carrasco et al 2022 (In prep.)
earlytransplant_startdates_table = '../rice-mapping-japan/data/earlytransplanting_startdates.csv'
earlytransplant_startdates = pd.read_csv(earlytransplant_startdates_table, sep=",")
earlytransplant_enddates_table = '../rice-mapping-japan/data/earlytransplanting_enddates.csv'
earlytransplant_enddates = pd.read_csv(earlytransplant_enddates_table, sep=",")

secondtransplant_startdates_table = '../rice-mapping-japan/data/secondtransplanting_startdates.csv'
secondtransplant_startdates = pd.read_csv(secondtransplant_startdates_table, sep=",")
secondtransplant_enddates_table = '../rice-mapping-japan/data/secondtransplanting_enddates.csv'
secondtransplant_enddates = pd.read_csv(secondtransplant_enddates_table, sep=",")

print(earlytransplant_startdates)
print(earlytransplant_enddates)
print(secondtransplant_startdates)
print(secondtransplant_enddates)<jupyter_output>Prefecture  8589  9094  9599  0004  0509  1014  1519
0      Kochi   102    97    95    94    94    94    93
1   Kumamoto   101    98   100   100   100   100   100
2   Miyazaki    91    85    82    80    80    80    79
3  Kagoshima    93    86    84    82    85    85    78
  Prefecture  8589  9094  9599  0004  0509  1014  1519
0      Kochi   120   119   117   114   112   112   112
1   Kumamoto   114   112   113   113   113   113   113
2   Miyazaki   105    99    97    95    95    97    96
3  Kagoshima   107   105   102   100   102   105   106
  Prefecture  8589  9094  9599  0004  0509  1014  1519
0      Kochi   218   218   218   218   218   218   218
1  Kagoshima   215   215   215   215   215   215   215
  Prefecture  8589  9094  9599  0004  0509  1014  1519
0      Kochi   228   228   228   228   228   228   228
1  Kagoshima   230   230   230   230   230   230   230<jupyter_text>Temperature data<jupyter_code># List of prefectures with DOY per temperature data
# These data can be found in Carrasco et al 2022 (In prep.)
temperature_table = '../rice-mapping-japan/data/prefectureDOYmintemp.csv'
tempdata = pd.read_csv(temperature_table, sep=",")
tempdata.head()<jupyter_output><empty_output><jupyter_text>Creates table with DOY with first/last 0/5/10 degrees<jupyter_code># Vectors to store TGS values (first day of certain min temperature)
TGS0S = []; TGS0E = []; TGS5S = []; TGS5E = []; TGS10S = []; TGS10E = [];


# We loop through each prefecture
for i in range(len(tempdata)):
    #print(i)
#i = 0

    # Extracts first DOY with mintemp > 0,5,10
    T0S = next(x for x, val in enumerate(tempdata.iloc[i,2:]) if val > 0)
    T5S = next(x for x, val in enumerate(tempdata.iloc[i,2:]) if val > 5)
    T10S = next(x for x, val in enumerate(tempdata.iloc[i,2:]) if val > 10)

    # Extracts first DOY with mintemp < 0,5,10 after summer (we consider temp values since August (DOY=214))
    # If value does not exist (desired temp falls into next year), we set to last day of year (so to not to use following year image)
    T0E = next((x for x, val in enumerate(tempdata.iloc[i,215:]) if val < 0), 365-214) + 214
    T5E = next((x for x, val in enumerate(tempdata.iloc[i,215:]) if val < 5), 365-214) + 214
    T10E = next((x for x, val in enumerate(tempdata.iloc[i,215:]) if val < 10), 365-214) + 214

    TGS0S.append(T0S); TGS5S.append(T5S); TGS10S.append(T10S); TGS0E.append(T0E); TGS5E.append(T5E); TGS10E.append(T10E); 

tgstable = pd.DataFrame({'Prefecture': tempdata['Prefecture'], 'TGS0S': TGS0S, 'TGS0E': TGS0E, 'TGS5S': TGS5S
                        , 'TGS5E': TGS5E, 'TGS10S': TGS10S, 'TGS10E': TGS10E})
tgstable.head()
# Kagoshima gives wrong values due to the TGS5S (above 5 in January but drops later), so we mannualy edit to real value
tgstable['TGS5S'][45] = 40
#tgstable.loc[tgstable['Prefecture']== 'Kagoshima']['TGS5S'] = 40
tgstable.loc[45]<jupyter_output><ipython-input-15-37567713f0f6>:2: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  tgstable['TGS5S'][45] = 40<jupyter_text>Load and process datasets for masks 3 and 4<jupyter_code># For mask 3: Forest (based on JAXA map)
# Load JAXA map
lc2015_path = 'users/luiscartor/jaxaLC2015'
lc2015 = geemap.ee.Image(lc2015_path)
# Select forest class
forest = lc2015.updateMask(lc2015.gt(5)).updateMask(lc2015.lt(10)).gt(5).selfMask()


# For mask 4: Elevation data
# Call SRTM elevation dataset
elevdataset = ee.Image('JAXA/ALOS/AW3D30/V2_2')
elevation = elevdataset.select('AVE_DSM')
# Calculate slope
slope = ee.Terrain.slope(elevation)
# Calculate slope > 20 deg
slopegt3 = slope.gt(20) 

# Uncomment to export masks
#geemap.ee_export_image_to_drive(forest, description='forestmask', folder='GEEexports', region=japan.geometry(), scale=30)<jupyter_output><empty_output><jupyter_text>2. MAIN ROUTINE  Extract collections<jupyter_code># Landsat 5
L5col = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR') \
    .select('B[1-7]','pixel_qa') \
    .map(addyear) \
    .map(addDOY)
       
# Landsat 7
L7col = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR') \
    .select('B[1-7]','pixel_qa') \
    .map(addyear) \
    .map(addDOY)

# Merge collections
L57col = L5col.merge(L7col)<jupyter_output><empty_output><jupyter_text>5 Year period collections<jupyter_code># Subset to 5-year periods
L57_8589 = L57col.filterDate('1985-01-01', '1989-12-31')
L57_9094 = L57col.filterDate('1990-01-01', '1994-12-31')
L57_9599 = L57col.filterDate('1995-01-01', '1999-12-31')
L57_0004 = L57col.filterDate('2000-01-01', '2004-12-31')
L57_0509 = L57col.filterDate('2005-01-01', '2009-12-31')
L57_1014 = L57col.filterDate('2010-01-01', '2014-12-31')
L57_1519 = L57col.filterDate('2015-01-01', '2019-12-31')

# Creates dictionary with periods and collections
periods_dict = {'8589': L57_8589,
              '9094': L57_9094,
              '9599': L57_9599,
              '0004': L57_0004,
              '0509': L57_0509,
              '1014': L57_1014,
              '1519': L57_1519}<jupyter_output><empty_output><jupyter_text>Main loop<jupyter_code># Creates empty list for storing national rice maps for each period
periodmapslist = list()

# LOOPS THROUGH EACH PERIOD
for periodcol in periods_dict.values():
# INDENT from here


# Uncomment and unindent to run for single period
#periodcol = L57_1519

    # Empty list to store prefectural maps for current loop period
    prefmapslist = list()


    # LOOPS THROUGH EACH PREFECTURE
    for prefecture in transplant_startdates['Prefecture']:
    # INDENT from here


        # Extracts prefecture polygon
        prefpolygon = prefs \
            .filter(ee.Filter.eq("NAME_1", prefecture)) 
        # Clips collection to prefecture
        prefcol = periodcol.filterBounds(prefpolygon)



        # PRE-PROCESSING
        
        # Cloud masking
        prefcol_masked = prefcol.map(landsatmasking)
        
        # Add Indices
        prefcol_indices = prefcol_masked.map(addINDICES)
        # Subset to spectral indices
        prefcol_ind = prefcol_indices.select('NDVI','LSWI','NDSI','EVI')



        # EXTRACTS FLOODING PERIODS
        
        # Flooding start day is considered the first day of transplanting (TS)
        flood_start = int(transplant_startdates[list(periods_dict.keys())[list(periods_dict.values()).index(periodcol)]][transplant_startdates[transplant_startdates['Prefecture']==prefecture].index[0]])
        # Flooding last day is considered the end day of transplanting + 30 days
        flood_end = int(transplant_enddates[list(periods_dict.keys())[list(periods_dict.values()).index(periodcol)]][transplant_enddates[transplant_enddates['Prefecture']==prefecture].index[0]] + 30)

        # Subset to flooding period
        prefcol_flood = prefcol_ind.filter(ee.Filter.gte('DOY', flood_start)).filter(ee.Filter.lte('DOY', flood_end))

        # Subset to non-flooding period
        prefcol_noflood1 = prefcol_ind.filter(ee.Filter.lte('DOY', flood_start))
        prefcol_noflood2 = prefcol_ind.filter(ee.Filter.gte('DOY', flood_end))
        prefcol_noflood = prefcol_noflood1.merge(prefcol_noflood2)        




        # EXTRACT FLOODING PERIODS FOR EARLY AND SECOND TRANSPLANTATIONS (ONLY USED FOR SOME PREFECTURES)
        
        # Prefectures with early transplantation periods
        if (prefecture =='Kochi' or prefecture =='Kumamoto' or prefecture =='Miyazaki' or prefecture =='Kagoshima'): 
            earlyflood_start = int(earlytransplant_startdates[list(periods_dict.keys())[list(periods_dict.values()).index(periodcol)]][earlytransplant_startdates[earlytransplant_startdates['Prefecture']==prefecture].index[0]])
            earlyflood_end = int(earlytransplant_enddates[list(periods_dict.keys())[list(periods_dict.values()).index(periodcol)]][earlytransplant_enddates[earlytransplant_enddates['Prefecture']==prefecture].index[0]] + 30)
            earlyprefcol_flood = prefcol_ind.filter(ee.Filter.gte('DOY', earlyflood_start)).filter(ee.Filter.lte('DOY', earlyflood_end))

        # Prefectures with second transplantation periods
        if (prefecture =='Kochi' or prefecture =='Kagoshima'): 
            secondflood_start = int(secondtransplant_startdates[list(periods_dict.keys())[list(periods_dict.values()).index(periodcol)]][secondtransplant_startdates[secondtransplant_startdates['Prefecture']==prefecture].index[0]])
            secondflood_end = int(secondtransplant_enddates[list(periods_dict.keys())[list(periods_dict.values()).index(periodcol)]][secondtransplant_enddates[secondtransplant_enddates['Prefecture']==prefecture].index[0]] + 30)
            secondprefcol_flood = prefcol_ind.filter(ee.Filter.gte('DOY', earlyflood_start)).filter(ee.Filter.lte('DOY', earlyflood_end))


        # FILTER COLLECTIONS FOR TEMPERATURE PERIODS
        
        # Extracting DOY for temperature data
        TGS0S = int(tgstable['TGS0S'][tgstable[tgstable['Prefecture']==prefecture].index[0]])
        TGS0E = int(tgstable['TGS0E'][tgstable[tgstable['Prefecture']==prefecture].index[0]]) 
        TGS5S = int(tgstable['TGS5S'][tgstable[tgstable['Prefecture']==prefecture].index[0]]) 
        TGS5E = int(tgstable['TGS5E'][tgstable[tgstable['Prefecture']==prefecture].index[0]]) 
        TGS10S = int(tgstable['TGS10S'][tgstable[tgstable['Prefecture']==prefecture].index[0]]) 
        TGS10E = int(tgstable['TGS10E'][tgstable[tgstable['Prefecture']==prefecture].index[0]])


        # Subset to period between first 0 and last 0 min degrees:
        L57_0S0E = prefcol_ind.filter(ee.Filter.gte('DOY', TGS0S)).filter(ee.Filter.lte('DOY', TGS0E))
        # Subset to period between first 5 and last 5 min degrees:
        L57_5S5E = prefcol_ind.filter(ee.Filter.gte('DOY', TGS5S)).filter(ee.Filter.lte('DOY', TGS5E))
        # Subset to period between first 0 and first 10 min degrees:
        L57_0S10S = prefcol_ind.filter(ee.Filter.gte('DOY', TGS0S)).filter(ee.Filter.lte('DOY', TGS10S))
        # Subset to period between first 10 + 40days and last 10 min degrees:
        L57_10S40to5E = prefcol_ind.filter(ee.Filter.gte('DOY', TGS10S+40)).filter(ee.Filter.lte('DOY', TGS10E))
        # Subset to period between last 10 and first 10 (not thermal growing season)
        L57_10E10S1 = prefcol_ind.filter(ee.Filter.lte('DOY', TGS10S))
        L57_10E10S2 = prefcol_ind.filter(ee.Filter.gte('DOY', TGS10E))
        L57_10E10S = L57_10E10S1.merge(L57_10E10S2)


        # Applies function to decide percentile based on transplantation period length
        percentile = percentileextract(flood_end-flood_start)
        L57_flood_lswi = prefcol_flood.select('LSWI').reduce(ee.Reducer.percentile([percentile])).rename('LSWI')
        L57_flood_ndvi = prefcol_flood.select('NDVI').median()
        L57_flood_evi = prefcol_flood.select('EVI').median()

        L57_flood = L57_flood_lswi.addBands(L57_flood_ndvi).addBands(L57_flood_evi)


        
        # APPLY MASKS (modified from on Dong 2016, RSE)

        # Mask 1) Sparce vegetation: soil, built-up, water body, low vegetated lands. 
        L57_5S5E_max = L57_5S5E.max()
        mask1 = L57_5S5E_max.select('EVI').lt(0.5)  

        # Mask 2) Natural vegetation mask: forests, natural wetlands, grass. 
        L57_10E10S_median = L57_10E10S.median()
        mask2 = L57_10E10S_median.select('EVI').gt(0.4) 

        # Mask 3) Forest mask 
        mask3 = forest

        # Mask 4) Slope
        mask4 = slopegt3

        

        # RICE CLASSIFICATION

        # Apply masks
        L57_flood_masked = landcovermasking(L57_flood)
        # Classify rice
        riceadded = addRICE(L57_flood_masked).select('RICE')
        riceclass = classornot(riceadded)



        # CLASSIFY RICE FOR PREFECTURES WITH EARLY AND SECOND PLANTATIONS
        # Prefectures with early transplantation periods
        if (prefecture =='Kochi' or prefecture =='Kumamoto' or prefecture =='Miyazaki' or prefecture =='Kagoshima'): 

            # Classify rice
            earlypercentile = percentileextract(earlyflood_end-earlyflood_start)
            earlyprefcol_flood_lswi = earlyprefcol_flood.select('LSWI').reduce(ee.Reducer.percentile([earlypercentile])).rename('LSWI')
            earlyprefcol_flood_ndvi = earlyprefcol_flood.select('NDVI').median()
            earlyprefcol_flood_evi = earlyprefcol_flood.select('EVI').median()
            earlyprefcol_flood = earlyprefcol_flood_lswi.addBands(earlyprefcol_flood_ndvi).addBands(earlyprefcol_flood_evi)

            # Apply mask
            earlyprefcol_flood_masked = landcovermasking(earlyprefcol_flood)

            # Add early rice to rice map
            earlyrice = addRICE(earlyprefcol_flood_masked).select('RICE') 
            earlyriceclass = classornot(earlyrice).multiply(10)
            riceclass = riceclass.add(earlyriceclass)




        # Prefectures with second transplantation periods
        if (prefecture =='Kochi' or prefecture =='Kagoshima'): 

            # Classify rice
            secondpercentile = percentileextract(secondflood_end-secondflood_start)
            secondprefcol_flood_lswi = secondprefcol_flood.select('LSWI').reduce(ee.Reducer.percentile([secondpercentile])).rename('LSWI')
            secondprefcol_flood_ndvi = secondprefcol_flood.select('NDVI').median()
            secondprefcol_flood_evi = secondprefcol_flood.select('EVI').median()
            secondprefcol_flood = secondprefcol_flood_lswi.addBands(secondprefcol_flood_ndvi).addBands(secondprefcol_flood_evi)

            # Apply masks
            secondprefcol_flood_masked = landcovermasking(secondprefcol_flood)

            # Add second rice to rice map
            secondrice = addRICE(secondprefcol_flood_masked).select('RICE') 
            secondriceclass = classornot(secondrice).multiply(20)
            riceclass = riceclass.add(secondriceclass)



        # CONVERT MAP TO BINARY
        # If parameter binarymap is 'yes', the map is converted into binary rice map
        if (binarymap == 'yes'):
            riceclass = riceclass.gte(0.1) 


        # MOSAICKING
        # Clip image to prefecture
        riceclass_clipped = riceclass.clip(prefpolygon)

        # Insert map in list
        prefmapslist.append(riceclass_clipped)


        # END OF PREFECTURES LOOP

    # Mosaic all prefectural maps for current loop period
    periodricecol = ee.ImageCollection.fromImages(prefmapslist)
    periodmosaic = periodricecol.mosaic()


    # Store period mosaic in list
    periodmapslist.append(periodmosaic)


    # FILTERING
    # It is possible to apply a spatial filtering (not used if not exported)
    periodmosaic_filtered = periodmosaic.reduceNeighborhood(**{   
      'reducer': ee.Reducer.mode(),   
      'kernel': ee.Kernel.square(1), 
      })


    # EXPORT MAPS
    # Exports to drive: edit name
    filename = 'tambo_'+list(periods_dict.keys())[list(periods_dict.values()).index(periodcol)]
    geemap.ee_export_image_to_drive(periodmosaic, description=filename, folder='GEEexports', region=japan.geometry(), scale=30)<jupyter_output>Exporting tambo_8589 ...
Exporting tambo_9094 ...
Exporting tambo_9599 ...
Exporting tambo_0004 ...
Exporting tambo_0509 ...
Exporting tambo_1014 ...
Exporting tambo_1519 ...
{'fraction_non_alphanumeric': 0.11656125556336379, 'fraction_numerical': 0.05598500819864137, 'mean_word_length': 3.901492537313433, 'pattern_counts': {'":': 0, '<': 56, '<?xml version=': 0, '>': 53, 'https://': 0, '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': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '130067', 'license': 'CC0-1.0', 'n_tokens_mistral': 8324, 'n_tokens_neox': 7340, 'n_words': 1586, 'path': 'ricemapping_main.ipynb', 'repo_name': 'luiscartor/rice-mapping-japan', 'size': '21345'}