<jupyter_start><jupyter_text>Creating random reference data polygonsUsing Collect Earth Online (CEO), these randomly scattered polygons will be classified as crop, non-crop, or mixed and form the knowledge base for the machine learning classifier.Stratifying points using the GFSAD cropland extent map, resampled to 60m resolution.        import os    os.chdir('data/GFSAD/')    gdalwarp -tr 0.0005389891704717132 0.0005389891704717132  -r mode GFSAD_mosaic.tif GFSAD_mosaic_60m.tif> TODO: Add an extra column identifying the GFSAD class (not just whether it came from the dataset). Use the label `GFSAD_CLASS`  Load modules<jupyter_code>import geopandas as gpd
from datacube.utils.cog import write_cog
import xarray as xr
import numpy as np
import pandas as pd
import datacube
from affine import Affine
from datacube.utils.geometry import assign_crs

import sys
sys.path.append('../Scripts')
from deafrica_spatialtools import xr_rasterize
from deafrica_dask import create_local_dask_cluster

import warnings
warnings.filterwarnings("ignore")
create_local_dask_cluster()<jupyter_output><empty_output><jupyter_text>Analysis Parameters<jupyter_code># asap_crop_mask = 'data/asap_mask_crop_v03.tif'
crop_mask_shp = 'data/GFSAD/GFSAD_mosaic_60m.tif'

#AEZ to mask crop mask
aez = 'data/AEZs/Indian_ocean.shp'

#location and file name to store shapefiles.
results = 'data/training_validation/collect_earth/indian_ocean/indian_ocean_TD'

#total number of sample points
n_sample = 1000

#number of classes to sample
n_class= 2<jupyter_output><empty_output><jupyter_text>Open datasets<jupyter_code>#open AEZ shapefile
aez = gpd.read_file(aez)

#find bouding box
xmin,ymin,xmax,ymax = aez.bounds.values[0][0], aez.bounds.values[0][1], aez.bounds.values[0][2], aez.bounds.values[0][3]

#plot to confirm region
aez.plot()
#open GFSAD cropmask at bounding coords of AEZ
da = xr.open_rasterio(crop_mask_shp).squeeze().sel(x=slice(xmin,xmax),
                                                  y=slice(ymax,ymin)).chunk({'x':3000, 'y':3000})
attrs = da.attrs

#rename x and y because for some reason the geobox params are different than xarray coords
da = da.rename({'x':'longitude', 'y':'latitude'})<jupyter_output><empty_output><jupyter_text>Mask by AEZ<jupyter_code>aez_raster = xr_rasterize(aez, da)
da = da.where(aez_raster).astype(np.int8)
da = da.to_dataset(name='GFSAD')
da.attrs = attrs
da = da.compute()<jupyter_output><empty_output><jupyter_text>Generate random points2000 random points: 50% (1000) in 2, 25% (500) in 1, 25% (500) in 02=crop1=non-crop0=no-data<jupyter_code>#rename back to x & y so code below works
da = da.rename({'longitude':'x', 'latitude':'y'})
class_sizes =[]
for class_id in np.arange(1, n_class+1):
    print(class_id)
    class_sizes.append((da.GFSAD==class_id).sum().values)

class_sizes = np.array(class_sizes)
print(class_sizes)
print(class_sizes/class_sizes.sum())
n_sample_class = np.ceil(n_sample*1./ n_class).astype(int)
print(n_sample_class)
label_picked = {}
for class_id in np.arange(1, n_class+1):
    if class_sizes[class_id-1]> 1e9:
        # slightly over sample
        n_sample_over = np.ceil(n_sample_class*len(da.x)*len(da.y)/class_sizes[class_id-1]).astype(int)
        random_x = np.random.choice(np.arange(len(da.x)), n_sample_over, replace=False)
        random_y = np.random.choice(np.arange(len(da.y)), n_sample_over, replace=False)
        match = da.GFSAD.values[random_y, random_x] == class_id
        random_y, random_x = random_y[match], random_x[match]
        if len(random_y) < n_sample_class:
            print("Not enough points are picked, try increase the number of random points")
            break
        else:
            pick = np.random.choice(np.arange(len(random_y)), n_sample_class, replace=False)
            y, x = random_y[pick], random_x[pick]
    else:
        index = np.argwhere(da.GFSAD.values.flatten() == class_id).squeeze()
        picked = np.random.choice(index, n_sample_class, replace=False)
        # convert back to x, y 
        y, x  = np.unravel_index(picked, da.GFSAD.values.shape)
    label_picked[class_id] = (y, x)
    np.savetxt(f'{results}_class_{class_id}.csv', 
               np.vstack((da.y[y].values, 
                          da.x[x].values)).transpose(),fmt='%d', delimiter=',')

for class_id in np.arange(1, n_class+1):
    y, x = label_picked[class_id]
    df = pd.DataFrame({'y': da.y[y].values, 'x':da.x[x].values})
    df['class']=class_id
    if class_id ==1: 
        dfs = df
    else: 
        dfs = dfs.append(df, ignore_index=True)
gdf_points = gpd.GeoDataFrame(
            dfs,
            crs=da.crs,
            geometry=gpd.points_from_xy(dfs.x, dfs.y)).reset_index()

gdf_points = gdf_points.drop(['x', 'y'],axis=1)

gdf_points.plot(figsize=(8,6), column='class', legend=True)<jupyter_output><empty_output><jupyter_text>Add subsample of points from GFSADThis will serve as an internal validation dataset.  Append 100 samples from GFSAD per region, 50 crop & 50 non-crop, and shuffle them within the larger dataset.We need a row identifying them so we can extract after the classification.<jupyter_code>#open data
gfsad_shp = 'data/training_validation/GFSAD2015/cropland_prelim_validation_GFSAD.shp'
gfsad_gdf = gpd.read_file(gfsad_shp)

#clip to AEZ
gfsad_gdf = gpd.overlay(gfsad_gdf,aez,how='intersection')
print('n crop: '+str(len(gfsad_gdf[gfsad_gdf['class'] ==1])))
gfsad_gdf.plot(column='class');
#sample 50 points from each class
crops = gfsad_gdf[gfsad_gdf['class']==1].sample(n=30, replace=False).reset_index(drop=True)
noncrops = gfsad_gdf[gfsad_gdf['class']==0].sample(n=30, replace=False).reset_index(drop=True)

#reclassify class int to match GFSAD scheme 
crops['class'] = np.where(crops['class'] == 1,2, crops['class'])
noncrops['class'] = np.where(noncrops['class'] == 0,1, crops['class'])
# add a column to identify validation samples
crops['GFSAD_sample']=True
noncrops['GFSAD_sample']=True
gdf_points['GFSAD_sample']=False
# Add GFSAD samples to dataset
gdf_points = pd.concat([gdf_points,crops,noncrops]).reset_index(drop=True).drop(columns=['ID','CODE','COUNTRY','index'])<jupyter_output><empty_output><jupyter_text>Create another small sample for training/testing<jupyter_code>crops = gfsad_gdf[gfsad_gdf['class']==1].sample(n=30, replace=False).reset_index(drop=True)
noncrops = gfsad_gdf[gfsad_gdf['class']==0].sample(n=30, replace=False).reset_index(drop=True)

#reclassify class int to match GFSAD scheme 
crops['class'] = np.where(crops['class'] == 1,2, crops['class'])
noncrops['class'] = np.where(noncrops['class'] == 0,1, crops['class'])

testing_gdf_points = pd.concat([crops,noncrops]).reset_index(drop=True).drop(columns=['ID','CODE','COUNTRY'])
testing_gdf_points.head()<jupyter_output><empty_output><jupyter_text>Shuffle rows so classes aren't sequential in CEO<jupyter_code>gdf_points = gdf_points.sample(frac=1).reset_index(drop=True)
testing_gdf_points = testing_gdf_points.sample(frac=1).reset_index(drop=True)<jupyter_output><empty_output><jupyter_text>Add a PLOTID field to satisfy Collect Earth<jupyter_code>gdf_points['PLOTID'] = range(0,len(gdf_points))
gdf_points['SAMPLEID'] = range(0,len(gdf_points))

testing_gdf_points['PLOTID'] = range(0,len(testing_gdf_points))
testing_gdf_points['SAMPLEID'] = range(0,len(testing_gdf_points))

gdf_points.head()<jupyter_output><empty_output><jupyter_text>Convert points to square polygons<jupyter_code>#set radius (in metres) around points
radius = 20

#convert to equal area to set polygon size in metres
gdf_poly = gdf_points.to_crs('EPSG:6933')
testing_gdf_poly = testing_gdf_points.to_crs('EPSG:6933')

#create circle buffer around points, then find envelope
gdf_poly['geometry'] = gdf_poly['geometry'].buffer(radius).envelope
testing_gdf_poly['geometry'] = testing_gdf_poly['geometry'].buffer(radius).envelope

#Convert back to lat/lon
gdf_poly = gdf_poly.to_crs('EPSG:4326')
testing_gdf_poly = testing_gdf_poly.to_crs('EPSG:4326')
#export to file
gdf_poly.to_file(results + '_polys.shp')
gdf_points.to_file(results + '_points.shp')

testing_gdf_poly.to_file(results + '_testing_polys.shp')
testing_gdf_points.to_file(results + '_testing_points.shp')<jupyter_output><empty_output><jupyter_text>Split datasets in half<jupyter_code>import geopandas as gpd
x=gpd.read_file('data/training_validation/collect_earth/northern/northern_TD_polys.shp')
print(len(x))
one = x[0:1100]
two = x[1100:]
one.to_file('data/training_validation/collect_earth/northern/northern_TD_polys_haithem.shp')
two.to_file('data/training_validation/collect_earth/northern/northern_TD_polys_malack.shp')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13160645835753282, 'fraction_numerical': 0.02067603670577303, 'mean_word_length': 5.844197138314786, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 47, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '378391', 'license': 'Apache-2.0', 'n_tokens_mistral': 3361, 'n_tokens_neox': 3215, 'n_words': 598, 'path': 'pre-post_processing/create_reference_data_polygons.ipynb', 'repo_name': 'eefaye/crop-mask', 'size': '8609'}