<jupyter_start><jupyter_text>Preparation of a labeled set of aerial images from public data sourcesThis notebook illustrates how a set of aerial images, labeled by land use classification, was generated from freely available U.S. datasets. The image sets generated in this notebook were used to train a DNN for image classification and evaluate its performance. For more detail, please see the rest of the [Embarrassingly Parallel Image Classification](https://github.com/Azure/Embarrassingly-Parallel-Image-Classification) repository.This walkthrough will demonstrate the technique using source data from Middlesex County, MA, the home of Microsoft's New England Research & Development (NERD) Center. The method was applied to many U.S. counties in parallel to generate the full training and validation sets.Note: it is not necessary to complete all steps in this tutorial in order to proceed to other tutorials in this series. For a fast start, complete the steps under the [Prepare an Azure Data Science Virtual Machine for image extraction](dsvm) section, then proceed directly to the [Dataset preparation for deep learning](prep) section.  Outline- [Source data](sourcedata)   - [National Agriculture Imagery Program](naip)   - [National Land Cover Database](nlcd)- [Prepare an Azure Data Science Virtual Machine for image extraction](dsvm)   - [Provision the VM](provision)   - [Connect to the VM by remote desktop](rd)   - [Install Python packages and supporting software](install)   - [Download input data locally](download)- [Produce image sets for training and evaluation](production)   - [Overview of the approach](overview)   - [Convert NAIP images from MrSID to GeoTIFF](conversion)   - [Add tiling for fast data extraction](tiling)   - [Find a lat-lon bounding box for the region](box)   - [Get the approximate dimensions of the bounding box in meters](meters)   - [Find tiles with consistent land use class](consistent)   - [Extract and save images for tiles of interest](extraction)   - [Dataset partitioning](partition)- [Dataset preparation for deep learning](prep)   - [Cognitive Toolkit (CNTK)](cntk)   - [TensorFlow](tf)- [Next Steps](nextsteps)  Source data National Agriculture Imagery ProgramThe US [National Agriculture Imagery Program](https://www.fsa.usda.gov/programs-and-services/aerial-photography/imagery-programs/naip-imagery/index), run by the [US Department of Agriculture](https://www.usda.gov/)'s [Farm Service Agency](https://www.fsa.usda.gov/), provides images within one month of collection every 2 years, in natural (RGB) color with 1-meter ground sample distance. In this tutorial, we use Compressed County Mosaics obtained from the [Geospatial Data Gateway](https://gdg.sc.egov.usda.gov/); you can download the 2016 Middlesex County, MA directly from [our mirror](https://mawahstorage.blob.core.windows.net/aerialimageclassification/naipsample/ortho_imagery_NAIPM16_ma017_3344056_01.zip) if you prefer.A zoomed-out view of the 2016 NAIP data covering Middlesex County, MA is shown below (Cambridge and Boston at ESE edge):  National Land Cover DatabaseThe US [National Land Cover Database](https://www.mrlc.gov/nlcd2011.php) (NLCD), maintained by the [Multi-Resolution Land Characteristic Consortium](https://www.mrlc.gov/), provides land use labels at 30-meter spatial resolution for the United States. The [sixteen land use classes](https://www.mrlc.gov/nlcd11_leg.php) in the NLCD are organized hierarchically into major (Developed, Forested, Cultivated, &c.) and minor (e.g. Deciduous, Evergreen, or Mixed Forest) categories. Classifications are based largely on seasonally-collected satellite imagery, including data collected outside the visual spectrum. Because of the extensive post-processing required to prepare the NLCD dataset, new versions are published approximately every five years, often with a multi-year delay between data collection and publication. For this project, we use the [NLCD 2011 Land Cover](https://www.mrlc.gov/nlcd11_data.php) dataset.The colorized NLCD data for Middlesex County, MA and surrounding areas are shown below. Developed land is shown in shades of red, water/wetlands in blue/cyan, forested land in green, and cultivated land in magenta.For more information on the NLCD, please see the following publication:<NAME>., <NAME>., <NAME>., <NAME>., <NAME>., <NAME>., <NAME>., <NAME>., <NAME>., and <NAME>., 2015, [Completion of the 2011 National Land Cover Database for the conterminous United States - Representing a decade of land cover change information](http://bit.ly/1K7WjO3). Photogrammetric Engineering and Remote Sensing, v. 81, no. 5, p. 345-354   Install Python packages and supporting software  Python packagesTo read the resulting GeoTIFF, we will use the GDAL package from within Python. As of this writing, Christoph Gohlke's [Unofficial Windows Binaries for Python Extension Packages](http://www.lfd.uci.edu/~gohlke/pythonlibs/) provide the easiest means to install the GDAL package, binaries, and library. We recommend that you download the Python 3.5 (`cp35`), `win_amd64` wheels for the following packages from Christoph's site:- [GDAL](http://www.lfd.uci.edu/~gohlke/pythonlibs/gdal)- [PyPROJ](http://www.lfd.uci.edu/~gohlke/pythonlibs/pyproj)- [Basemap](http://www.lfd.uci.edu/~gohlke/pythonlibs/basemap)To install these wheels:1. Launch an Anaconda Prompt.1. Type the following commands (in this order):  ```activate py35  pip install [GDAL wheel filepath]  pip install [PyPROJ wheel filepath]  pip install [Basemap wheel filepath]  ```      If you encounter a "BadZipFile" error when installing the Basemap wheel, you may need to upgrade `pip` first:  ```  python -m pip install --upgrade pip  ```  Produce image sets for training and evaluation  Overview of the approachOur goal is to extract many smaller aerial images ("tiles") from the NAIP image of the entire county. Each tile will be labeled with a land use classification determined from the NLCD data. Before launching into the technical details of this process, we provide an overview of the steps involved.  Choosing tile sizeResNet DNNs are traditionally trained with 224 px x 224 px images: this sets our minimum tile size at 224 meters x 224 meters (because our aerial imagery has approx. 1 meter resolution). A sample tile with these dimensions is shown below at full size.  Proposing candidate tilesWe divided the Middlesex County, MA region along a rectangular grid of tiles, as illustrated in the approximately 1 km x 1 km region shown below: To create the grid, we will need to find the latitude/longitude boundaries and dimensions in meters of the county-level NAIP image. We do this by converting the image to GeoTIFF format and using the Python GDAL package to extract the necessary information from the header. Once we have defined the grid, we can use GDAL to extract the image for each tile from the county-level image.  Assigning land use classes to tilesThe tiles described above often contain multiple types of land, making it difficult to assign a single land use label. To illustrate this problem, we overlaid the NLCD land use labels and the NAIP aerial image in a region centered on the [Boston Common](https://binged.it/2kTjwOe) (NAIP-only image shown above): A straightforward approach to avoid ambiguity when assigning labels would be to remove any tiles with multiple land uses. Unfortunately, we found this requirement too restrictive: the number of tiles retained by this method was too small for model training. We developed an arbitrary, more lenient criterion for labeling in which we sampled the land use labels at nine regularly-spaced points within each tile, as depicted below:If all nine points had the same land use class, we assigned that land use class as the tile's label. Images that could not be assigned a label by this method were not included in the training and validation sets.  Add tiling for fast image extractionTo speed up extraction of images from arbitrary regions in the 45 GB file, we add tiling to the image. Since we must write the tiled image before deleting the untiled version, we briefly double our disk space usage at this step.<jupyter_code>##### NOTE: NEED TO UPDATE PATHNAMES FOR YOUR SYSTEM
naip_dir = "C:\\Users\\ME\\Dropbox\\HCSproject\\data\\PoC\\app_kalbar_cntk\\"
import os
import subprocess
import gdal
filename_base = "app_kalbar_input_ndvi_s2"
nlcd_filepath = naip_dir + 'app_kalbar_remap.tif'
output_dir = naip_dir + "tiles\\"

img = gdal.Open(os.path.join(naip_dir, '{}.tif'.format(filename_base)))
out = os.path.join(naip_dir, '{}_tiled.tif'.format(filename_base))
translateoptions = gdal.TranslateOptions(gdal.ParseCommandLine("-co COMPRESS=LZW, TILED=YES"))
gdal.Translate(out, img, options = translateoptions)
naip_filepath = out<jupyter_output><empty_output><jupyter_text>Find a lat-lon bounding box for the regionWe now extract the county's boundaries in latitude/longitude coordinates from the GeoTIFF file. Note that because of the county's irregular shape, some regions inside this lat/lon bounding box will not contain any data. The bounding box is nonetheless useful for establishing the dimensions of a tiling grid.<jupyter_code>from osgeo import gdal
from gdalconst import *
import osr
from mpl_toolkits.basemap import Basemap
from collections import namedtuple

LatLonBounds = namedtuple('LatLonBounds', ['llcrnrlat', 'llcrnrlon', 'urcrnrlat', 'urcrnrlon'])

def get_bounding_box(naip_filepath):
    ''' Finds a bounding box for the NAIP GeoTIFF in lat/lon '''
    naip_image = gdal.Open(naip_filepath, GA_ReadOnly)
    naip_proj = osr.SpatialReference()
    naip_proj.ImportFromWkt(naip_image.GetProjection())
    naip_ulcrnrx, naip_xstep, _, naip_ulcrnry, _, naip_ystep = naip_image.GetGeoTransform()

    world_map = Basemap(lat_0=0,
                        lon_0=0,
                        llcrnrlat=-90, urcrnrlat=90,
                        llcrnrlon=-180, urcrnrlon=180,
                        resolution='c', projection='stere')
    world_proj = osr.SpatialReference()
    world_proj.ImportFromProj4(world_map.proj4string)
    ct_to_world = osr.CoordinateTransformation(naip_proj, world_proj)
    
    lats = []
    lons = []
    for corner_x, corner_y in [(naip_ulcrnrx, naip_ulcrnry),
                               (naip_ulcrnrx, naip_ulcrnry + naip_image.RasterYSize * naip_ystep),
                               (naip_ulcrnrx + naip_image.RasterXSize * naip_xstep,
                                naip_ulcrnry + naip_image.RasterYSize * naip_ystep),
                               (naip_ulcrnrx + naip_image.RasterXSize * naip_xstep, naip_ulcrnry)]:
        xpos, ypos, _ = ct_to_world.TransformPoint(corner_x, corner_y)
        lon, lat = world_map(xpos, ypos, inverse=True)
        lats.append(lat)
        lons.append(lon)

    return(LatLonBounds(llcrnrlat=min(lats),
                        llcrnrlon=min(lons),
                        urcrnrlat=max(lats),
                        urcrnrlon=max(lons)))

region_bounds = get_bounding_box(naip_filepath)
os.environ['PROJ_LIB']<jupyter_output><empty_output><jupyter_text>Get the approximate dimensions of the bounding box in metersHere we approximate the width/height of the bounding box in meters, using the fact that latitude does not change substantially on the county scale. This dimensions in meters are needed later to define a grid of 224 meter x 224 meter tiles.<jupyter_code>RegionSize = namedtuple('RegionSize', ['width', 'height'])  # in meters!
import numpy as np

def get_approx_region_size(region_bounds):
    ''' Returns the region width (at mid-lat) and height in meters'''
    mid_lat_radians = (region_bounds.llcrnrlat + region_bounds.urcrnrlat) * \
                      (np.pi / 360)
    earth_circumference = 6.371E6 * 2 * np.pi # in meters
    region_middle_width_meters = (region_bounds.urcrnrlon - region_bounds.llcrnrlon) * \
                                 earth_circumference * np.cos(mid_lat_radians) / 360
    region_height_meters = (region_bounds.urcrnrlat - region_bounds.llcrnrlat) * \
                           earth_circumference / 360
    return(RegionSize(region_middle_width_meters, region_height_meters))

approx_region_size = get_approx_region_size(region_bounds)
26700/30<jupyter_output><empty_output><jupyter_text>Find tiles with consistent land use class Define helper functionsThe function below defines helper functions that help us access data from both input images (the NAIP aerial image and the NLCD land use label image) using latitude and longitude coordinate systems.<jupyter_code>from PIL import Image

def create_helper_functions(region_bounds, nlcd_filepath, naip_filepath):
    ''' Makes helper functions to label points (NLCD) and extract tiles (NAIP) '''
    nlcd_image = gdal.Open(nlcd_filepath, GA_ReadOnly)
    nlcd_proj = osr.SpatialReference()
    nlcd_proj.ImportFromWkt(nlcd_image.GetProjection())
    nlcd_ulcrnrx, nlcd_xstep, _, nlcd_ulcrnry, _, nlcd_ystep = nlcd_image.GetGeoTransform()
    
    naip_image = gdal.Open(naip_filepath, GA_ReadOnly)
    naip_proj = osr.SpatialReference()
    naip_proj.ImportFromWkt(naip_image.GetProjection())
    naip_ulcrnrx, naip_xstep, _, naip_ulcrnry, _, naip_ystep = naip_image.GetGeoTransform()
    
    region_map = Basemap(lat_0 = (region_bounds.llcrnrlat + region_bounds.urcrnrlat)/2,
                         lon_0 = (region_bounds.llcrnrlon + region_bounds.urcrnrlon)/2,
                         llcrnrlat=region_bounds.llcrnrlat,
                         llcrnrlon=region_bounds.llcrnrlon,
                         urcrnrlat=region_bounds.urcrnrlat,
                         urcrnrlon=region_bounds.urcrnrlon,
                         resolution='c',
                         projection='stere')
    
    region_proj = osr.SpatialReference()
    region_proj.ImportFromProj4(region_map.proj4string)
    ct_to_nlcd = osr.CoordinateTransformation(region_proj, nlcd_proj)
    ct_to_naip = osr.CoordinateTransformation(region_proj, naip_proj)

    def get_nlcd_label(point):
        ''' Project lat/lon point to NLCD GeoTIFF; return label of that point '''
        basemap_coords = region_map(point.lon, point.lat)  # NB unusual argument order#       
        #x, y, _ = [int(i) for i in ct_to_nlcd.TransformPoint(*basemap_coords)]
####### NOTE: ORIGINAL CODE USED LINE ABOVE. ROUNDING SEEMED TO CAUSE PROBLEMS, BUT NOW LATER CODE WONT WORK. MAYBE ITS AN ISSUE
####### WITH THE SHIFT IN RESOLUTIONS FROM THE NAIP/NLDC APPLICATION TO OUR APPLICATION

        x, y, _ = [i for i in ct_to_nlcd.TransformPoint(*basemap_coords)]
        xoff = int(round((x - nlcd_ulcrnrx) / nlcd_xstep))
        yoff = int(round((y - nlcd_ulcrnry) / nlcd_ystep))
        label = int(nlcd_image.ReadAsArray(xoff=xoff, yoff=yoff, xsize=1, ysize=1))
        return(label)
        
    def get_naip_tile(tile_bounds, tile_size):
        ''' Check that tile lies within county bounds; if so, extract its image '''
        
        # Transform tile bounds in lat/lon to NAIP projection coordinates
        xmax, ymax = region_map(tile_bounds.urcrnrlon, tile_bounds.urcrnrlat)
        xmin, ymin = region_map(tile_bounds.llcrnrlon, tile_bounds.llcrnrlat)
        xstep = (xmax - xmin) / tile_size.width
        ystep = (ymax - ymin) / tile_size.height

        grid = np.mgrid[xmin:xmax:tile_size.width * 1j, ymin:ymax:tile_size.height * 1j]
        shape = grid[0, :, :].shape
        size = grid[0, :, :].size
        xy_target = np.array(ct_to_naip.TransformPoints(grid.reshape(2, size).T))
        xx = xy_target[:,0].reshape(shape)
        yy = xy_target[:,1].reshape(shape)
        
        # Extract rectangle from NAIP GeoTIFF containing superset of needed points
        xoff = int(round((xx.min() - naip_ulcrnrx) / naip_xstep))
        yoff = int(round((yy.max() - naip_ulcrnry) / naip_ystep))
        xsize_to_use = int(np.ceil((xx.max() - xx.min())/np.abs(naip_xstep))) + 1
        ysize_to_use = int(np.ceil((yy.max() - yy.min())/np.abs(naip_ystep))) + 1
        data = naip_image.ReadAsArray(xoff=xoff,
                                      yoff=yoff,
                                      xsize=xsize_to_use,
                                      ysize=ysize_to_use)        
        # Map the pixels of interest in NAIP GeoTIFF to the tile (might involve rotation or scaling)
        #image = np.zeros((xx.shape[1], xx.shape[0], 3)).astype(int)  # rows are height, cols are width, third dim is color
        image = np.zeros((xx.shape[1], xx.shape[0],1)).astype(int)  # rows are height, cols are width, third dim is color
        
        try:
            for i in range(xx.shape[0]):
                for j in range(xx.shape[1]):
                    x_idx = int(round((xx[i,j] - naip_ulcrnrx) / naip_xstep)) - xoff
                    y_idx = int(round((yy[i,j] - naip_ulcrnry) / naip_ystep)) - yoff
                    #image[xx.shape[1] - j - 1, i, :] = data[:, y_idx, x_idx]
                    image[xx.shape[1] - j - 1, i,:] = data[y_idx, x_idx]
        except TypeError as e:
            # The following can occur if our pixel superset request exceeds the GeoTIFF's bounds
            print("Out of Bounds")
            return(None)
        
        #if np.sum(image.sum(axis=2) == 0) > 10: # too many nodata pixels
        if np.sum(image.sum() == 0) > 10: # too many nodata pixels
            print("Too Many NoData")
            return None
        
        #image = Image.fromarray(image.astype('uint8'))
        image = Image.fromarray((image * 255).astype(np.uint8))
        return(image)
    
    return(get_nlcd_label, get_naip_tile)
get_nlcd_label, get_naip_tile = create_helper_functions(region_bounds, nlcd_filepath, naip_filepath)<jupyter_output><empty_output><jupyter_text>Tile selectionWe form a rectangular grid of candidate tiles spanning as much of the bounding box as possible. The boundaries of tiles that could be assigned a label are noted, and will later be used to extract the corresponding aerial image.<jupyter_code>LatLonPosition = namedtuple('LatLonPosition', ['lat', 'lon'])
Tile = namedtuple('Tile', ['bounds', 'label'])

nlcd_label_to_class = {0: 'NA',
                       1: 'Not_HCSA',
                       2: 'HCSA'}

def find_tiles_with_consistent_labels(region_bounds, region_size, tile_size):
    ''' Find tiles for which nine grid points all have the same label '''
    tiles_wide = int(np.floor(region_size.width / tile_size.width))
    tiles_tall = int(np.floor(region_size.height / tile_size.height))
    tile_width = (region_bounds.urcrnrlon - region_bounds.llcrnrlon) / (region_size.width / tile_size.width)
    tile_height = (region_bounds.urcrnrlat - region_bounds.llcrnrlat) / (region_size.height / tile_size.height)
    
    current_lat = region_bounds.llcrnrlat
    current_lon = region_bounds.llcrnrlon
    
    tiles_to_use = []
    for i in range(tiles_tall):
        for j in range(tiles_wide):
            try:
                labels = []
                for k in range(3):
                    for ell in range(3):
                        my_label = get_nlcd_label(LatLonPosition(lat=current_lat + tile_width * (1 + 2*k) / 6,
                                                                 lon=current_lon + tile_height * (1 + 2*ell) / 6))
                        labels.append(nlcd_label_to_class[my_label])
                num_matching = np.sum(np.array(labels) == labels[4])
                if (num_matching == 9):
                    bounds = LatLonBounds(llcrnrlat=current_lat,
                                          llcrnrlon=current_lon,
                                          urcrnrlat=current_lat + tile_height,
                                          urcrnrlon=current_lon + tile_width)
                    tiles_to_use.append(Tile(bounds=bounds,
                                             label=labels[4]))
            except KeyError:
                pass
            current_lon += tile_width
        current_lon = region_bounds.llcrnrlon
        current_lat += tile_height
    return(tiles_to_use)

tile_size = RegionSize(900, 900)
tiles = find_tiles_with_consistent_labels(region_bounds, approx_region_size, tile_size)
print('Found {} tiles to extract'.format(len(tiles)))<jupyter_output><empty_output><jupyter_text>Extract and save images for tiles of interestThe most time-consuming step is the extraction of images from the GeoTIFF. We use a helper function defined earlier, `naip_get_tile`, to check that the tile lies entirely inside the county shapefile and, if so, return the extracted image. Images are sorted into directories based on land use class. Expect this step to take several hours. You can monitor the progress by checking the number of images created in the output directory.<jupyter_code>def extract_tiles(tiles, dest_folder, filename_base):
    ''' Coordinates saving tile data, including extracted images and CSV descriptions '''
    if not os.path.exists(dest_folder):
        os.makedirs(dest_folder)
        
    tile_descriptions = []
    i = 0
    while i < len(tiles):
        tile = tiles[i]
        print(tile)
        tile_image = get_naip_tile(tile.bounds, tile_size)
        i += 1
        if (tile_image is None):
            continue  # tile did not lie entirely within the county boundary (it was at least partially blank)
        my_directory = os.path.join(dest_folder, '{}'.format(tile.label))
        my_filename = os.path.join(my_directory, '{}_{}.png'.format(filename_base, i))
        print(my_filename)
        if not os.path.exists(my_directory):
            os.makedirs(my_directory)
        tile_image.save(my_filename, 'PNG')
    return

extract_tiles(tiles, output_dir, filename_base)
filename_base


len(tiles)
        
    tile_descriptions = []
    i = 0
    while i < len(tiles):
        tile = tiles[i]
        print(tile)
        tile_image = get_naip_tile(tile.bounds, tile_size)
        i += 1
        if (tile_image is None):
            continue  # tile did not lie entirely within the county boundary (it was at least partially blank)
        my_directory = os.path.join(dest_folder, '{}'.format(tile.label))
        my_filename = os.path.join(my_directory, '{}_{}.png'.format(filename_base, i))
        if not os.path.exists(my_directory):
            os.makedirs(my_directory)
        tile_image.save(my_filename, 'PNG')<jupyter_output><empty_output><jupyter_text>Dataset partitioningTo create large and diverse training/validation sets, we repeated the steps above for 12 counties spread across the United States. We partitioned the images into training and validation sets at the county level: this division allows a reasonably accurate estimate of model performance by ensuring some dissimilarity between images seen during training and those used to evaluate the model's performance. We balanced each image set by randomly removing images from overrepresented classes.While you are welcome to repeat those steps yourself, for expediency we recommend that you continue with this tutorial by downloading our training and validation sets (linked in the next section).  Prepare deep learning framework-specific input filesIf you have not generated your own training and validation sets through image extraction, download the following files and decompress them in your VM's `D` drive:- [Balanced training image set (~3 GB)](https://mawahstorage.blob.core.windows.net/aerialimageclassification/imagesets/balanced_training_set.zip)- [Balanced validation image set (~1 GB)](https://mawahstorage.blob.core.windows.net/aerialimageclassification/imagesets/balanced_validation_set.zip)The image sets linked above contain raw PNG images sorted into folders by their assigned label. Many deep learning frameworks require proprietary image formats or supporting files to efficiently load images in minibatches for training. We now produce the supporting files needed by our CNTK and TensorFlow training scripts. Note that we do not need to produce similar supporting files for the validation image set, because we will not use minibatching when applying the trained models to the validation set.Update the `training_image_dir` variable below to reflect the directory where your training and validation sets have been saved. The `label_to_number_dict` variable specifies the correspondence between the label names and a numeric code; it does not need to be modified unless you have changed the labeling scheme.<jupyter_code>import os
import numpy as np
import pandas as pd
import tensorflow as tf

training_image_dir = 'C:\\Users\\ME\\Dropbox\\HCSproject\\data\\PoC\\app_kalbar_cntk\\tiles\\balanced_training_set'
label_to_number_dict = {'NA': 0,
                        'Not_HCSA': 1,
                        'HCSA': 2}<jupyter_output><empty_output><jupyter_text>Cognitive Toolkit (CNTK)Our CNTK training script uses a MAP file -- a tab-delimited file where each line lists an image's filepath and label -- to load image data in minibatches. We generate the MAP file as follows:<jupyter_code>with open(os.path.join(training_image_dir, 'map.txt'), 'w') as map_file:
    for label in np.sort(os.listdir(training_image_dir)):
        my_dir = os.path.join(training_image_dir, label)
        if not os.path.isdir(my_dir):
            continue
        for filename in os.listdir(my_dir):
            map_file.write('{}\t{}\n'.format(os.path.join(my_dir, filename),
                                             label_to_number_dict[label]))<jupyter_output><empty_output><jupyter_text>TensorFlowWe made use of the [`tf-slim` API](https://github.com/tensorflow/tensorflow/tree/master/tensorflow/contrib/slim) for TensorFlow, which provides pre-trained ResNet models and helpful scripts for retraining and scoring. These scripts require converting the image set into [TFRecords](https://www.tensorflow.org/how_tos/reading_data/file_formats) for minibatching. (Each TFRecord contains many image files as well as their labels.) We also create a `labels.txt` file mapping the labels to integer values, and a `dataset_split_info.csv` file describing the images assigned to the training set.The following code was modified from the [TensorFlow models repo's slim subdirectory](https://github.com/tensorflow/models/tree/master/slim).<jupyter_code># Original Copyright 2016 The TensorFlow Authors. All Rights Reserved.
# Modified 2017 by Microsoft Corporation.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ==============================================================================
np.random.seed(5318)

class ImageReader(object):
    def __init__(self):
        # Initializes function that decodes RGB JPEG data.
        self._decode_png_data = tf.placeholder(dtype=tf.string)
        self._decode_png = tf.image.decode_png(self._decode_png_data, channels=3)

    def read_image_dims(self, sess, image_data):
        image = self.decode_png(sess, image_data)
        return image.shape[0], image.shape[1]

    def decode_png(self, sess, image_data):
        image = sess.run(self._decode_png,
                         feed_dict={self._decode_png_data: image_data})
        assert len(image.shape) == 3
        assert image.shape[2] == 3
        return image
    
def image_to_tfexample(image_data, image_format, height, width, class_id):
    return tf.train.Example(features=tf.train.Features(feature={
        'image/encoded': tf.train.Feature(bytes_list=tf.train.BytesList(value=[image_data])),
        'image/format': tf.train.Feature(bytes_list=tf.train.BytesList(value=[image_format])),
        'image/class/label': tf.train.Feature(int64_list=tf.train.Int64List(value=[class_id])),
        'image/height': tf.train.Feature(int64_list=tf.train.Int64List(value=[height])),
        'image/width': tf.train.Feature(int64_list=tf.train.Int64List(value=[width])),
    }))

def find_images(image_dir):
    training_filenames = []
    for folder in os.listdir(image_dir):
        folder_path = os.path.join(image_dir, folder)
        if not os.path.isdir(folder_path):
            continue
        ''' This is a new directory/label -- consider all images inside it '''
        my_filenames = []
        for filename in os.listdir(folder_path):
            my_filenames.append(os.path.join(folder_path, filename))
        training_filenames.extend(my_filenames)
    print('Found {} training images'.format(len(training_filenames)))
    return(training_filenames)
    
def write_dataset(dataset_name, split_name, my_filenames,  image_dir, n_shards=5):
    num_per_shard = int(np.ceil(len(my_filenames) / n_shards))
    records = []
    with tf.Graph().as_default():
        image_reader = ImageReader()
        with tf.Session('') as sess:
            for shard_idx in range(n_shards):
                shard_filename = os.path.join(image_dir,
                                              '{}_{}_{:05d}-of-{:05d}.tfrecord'.format(dataset_name,
                                                                                       split_name,
                                                                                       shard_idx+1,
                                                                                       n_shards))
                with tf.python_io.TFRecordWriter(shard_filename) as tfrecord_writer:
                    for image_idx in range(num_per_shard * shard_idx,
                                           min(num_per_shard * (shard_idx+1), len(my_filenames))):
                        with open(my_filenames[image_idx], 'rb') as f:
                            image_data = f.read()
                        height, width = image_reader.read_image_dims(sess, image_data)
                        class_name = os.path.basename(os.path.dirname(my_filenames[image_idx]))
                        class_id = label_to_number_dict[class_name]
                        example = image_to_tfexample(image_data, b'png', height, width, class_id)
                        tfrecord_writer.write(example.SerializeToString())
                        records.append([dataset_name, split_name, my_filenames[image_idx], shard_idx,
                                        image_idx, class_name, class_id])
    df = pd.DataFrame(records, columns=['dataset_name', 'split_name', 'filename', 'shard_idx', 'image_idx',
                                        'class_name', 'class_id'])
    return(df)
 
training_image_dir = 'C:\\Users\\ME\\Dropbox\\HCSproject\\data\\PoC\\app_kalbar_cntk\\tiles\\balanced_training_set'
training_filenames = find_images(training_image_dir)
training_filenames = np.random.permutation(training_filenames)
df = write_dataset('aerial', 'train', training_filenames, training_image_dir, n_shards=50)
df.to_csv(os.path.join(training_image_dir, 'dataset_split_info.csv'), index=False)

with open(os.path.join(training_image_dir, 'labels.txt'), 'w') as f:
    for key, value in label_to_number_dict.items():
        f.write('{}:{}\n'.format(key, value))<jupyter_output>Found 3691 training images
{'fraction_non_alphanumeric': 0.0890397857689365, 'fraction_numerical': 0.007969905636317266, 'mean_word_length': 3.2719596895002043, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 48, 'https://': 16, 'lorem ipsum': 0, 'www.': 13, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 43, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '46716', 'license': 'Apache-2.0', 'n_tokens_mistral': 9374, 'n_tokens_neox': 8945, 'n_words': 2852, 'path': 'image_set_preparation_hcsa.ipynb', 'repo_name': 'eggenm/Embarrassingly-Parallel-Image-Classification', 'size': '31404'}