<commit_before>""" Extract Patches

This scripts reads in a shapefile and a raster image extracts patches to be used for machine learning algorithms.
If the shapefile consists of points, then a rectangular patch centered at each point is extracted. If the shapefile
consists of lines, then the first and second point are used to determine an orientation, and a patch is extracted
centered at the midpoint of the line and rotated so that the line would be horizontal, pointing in the +x direction.

The ouput folder will always contain a CSV file with the shapefile and raster name, as well as index of the feature in
the shapefile, the latitude and longitude, and image row and column of each feature.


Usage:
    extract_patches.py (-i <shapefile>)...  -r <raster> [--odir <outdir>] [options]
    extract_patches.py -h

Options:
    -i <shapefile>           A shapefile with features indicated. Wildcards are allowed.
    -r <raster>              A  raster source image to extract features from. No wildcards -- if you have multiple
                             raster files consider using `gdalbuildvrt` first to create a virtual mosaic.
    --odir <outdir>          The output folder. [default: ./patches/]
    --otif                   Outout TIF files [default]
    --ojpg                   Output JPG files
    --size <width>,<height>  The size of each extracted patch, in pixels. Please make sure there are no space characters
                             between the arguments unless you put it in quotes! [default: 64,64]
    --scale <scale>          The amount to scale each image by before extracting the patches (e.g.  if the raster is
                             in a higher resolution you may choose to scale it down). [default: 1.0]
    --csv <csvfile>          Create a CSV file with the center lat/lon, center pixel xy, angle, and the patch filename
                             as well as the name of the shapefile source and the raster spource image.
    --silent                 Suppress progress updates.
    --vclip <min,max>        Clip input values to `min` and `max`. Mainly useful when floating point output is not an
                             option (e.g. JPEG output).
    --vstretch <min,max>     Stretch output values to 'min' and 'max', after clipping. Mainly useful when floating point
                             output is not possible. Note that JPEG output must be in 0 to 1.

"""
import hashlib
import os

import numpy
import osgeo
import rasterio
from affine import Affine
from math import degrees, atan2, hypot, floor, ceil
from osgeo import ogr, gdalnumeric
from rasterio._io import RasterReader, RasterUpdater
#from skimage.transform import rescale
from scipy.ndimage import rotate, zoom
from skimage.util.dtype import img_as_ubyte

BREAKING_RELEASE = 0  # Big, rare, when you actually remove deprecated code.
NONBREAKING_RELEASE = 0  # Upload new binaries when you increment this, much less frequently then minor
MINOR_UPDATE = 1  # Increment for pullers
EXPERIMENTAL = 1

VERSION = '{}.{}.{}'.format(BREAKING_RELEASE, NONBREAKING_RELEASE, MINOR_UPDATE)
if EXPERIMENTAL:
    VERSION += 'a{}'.format(EXPERIMENTAL)
__version__ = VERSION

from glob import glob
from docopt import docopt
import logging
from progressbar import Percentage, Bar, ETA, ProgressBar


def collect_filenames(wildcard_list):
    logger = logging.getLogger(__name__)
    results = []
    for wildcard in wildcard_list:
        matching_files = glob(wildcard)
        if len(matching_files) is 0:
            logger.warning("No files matching input specification '{}'".format(wildcard))
        for filename in matching_files:
            results.append(filename)
    return results


def main():
    args = docopt(__doc__, version=VERSION)
    logging.basicConfig(level=logging.INFO)
    logger = logging.getLogger('extract_patches')
    logger.debug('input \n {}'.format(args))

    assert isinstance(logger, logging.Logger)

    shapefiles = collect_filenames(args['-i'])
    if len(shapefiles) == 0:
        logger.error('No matching shapefiles for inoput `{}`'.format(args['-i']))
        return

    raster = args['-r']

    try:
        size = [int(x) for x in args['--size'].split(',')]
        patch_width, patch_height = size
    except:
        logger.error("Unable to parse option '--size'")
        return

    try:
        scale = float(args['--scale'])
        assert scale > 0
    except:
        logger.error("Unable to parse option '--scale'")
        return

    silent = args['--silent']

    output_folder = args['--odir']
    try:
        if not os.path.isdir(output_folder):
            os.makedirs(output_folder)
    except:
        logger.error("Unable to find or create output directory `{}`".format(output_folder))
        return

    if args['--ojpg']:
        fmt = '.jpg'
    else:  # args['--otif']  (default)
        fmt = '.tif'

    clip  = args['--vclip'] is not None
    if clip:
        clipmin, clipmax = [float(x) for x in args['--vclip'].split(',')]

    stretch = args['--vstretch'] is not None
    if stretch:
        stretchmin, stretchmax = [float(x) for x in args['--vstretch'].split(',')]


    # Estimate number of shape features
    count = 0
    if not silent:
        pbar = ProgressBar(len(shapefiles), ['Counting Features:', Percentage(), ' ', Bar(), ' ', ETA()])
        pbar.start()
    for i, s in enumerate(shapefiles):
        vector = ogr.Open(s)
        layer = vector.GetLayer()
        count += layer.GetFeatureCount()
        if not silent:
            pbar.update(i)
    if not silent:
        pbar.finish()

    with rasterio.open(raster) as rf:
        assert isinstance(rf, RasterReader)
        affine = rf.affine
        geo_to_pixels = ~affine

        if not silent:
            pbar = ProgressBar(count, ['Exporting Patches:', Percentage(), ' ', Bar(), ' ', ETA()])
            pbar.start()
        for sf in shapefiles:
            vector = ogr.Open(sf)
            layer = vector.GetLayer()
            for f in layer:
                if not silent:
                    pbar.update(pbar.currval+1)
                geom = f.GetGeometryRef()
                points = geom.GetPoints()
                source = points[0]
                target = points[-1]
                sx, sy = geo_to_pixels * source
                tx, ty = geo_to_pixels * target
                if len(points) == 2:
                    cx, cy = (sx + tx) / 2, (sy + ty) / 2
                else:
                    cx, cy = points[1]
                dx, dy = (tx - sx), (ty - sy)
                theta = degrees(atan2(dy, dx))  # In PIXELS, CCW from +x. Not necessarily CCW from E (or CW from N)
                r1 = hypot(tx - cx, ty - cy)
                r2 = hypot(cx - sx, cy - sy)
                r1, r2 = max(r1, r2), min(r1, r2)  # For 3 points, we assume two radii. Else these are duplicates.
                gx, gy = affine * (cx, cy)  # Geographic coordinates (e.g. lat lon) of the center.

                # We read a square slightly larger than the scaled version of our patch, so that
                # we can safely rotate the raster without missing pixels in the corners.

                box_radius = hypot(patch_width, patch_height) / (2.0 * scale)
                x0, x1 = int(floor(cx - box_radius)), int(ceil(cx + box_radius))
                y0, y1 = int(floor(cy - box_radius)), int(ceil(cy + box_radius))

                                # save patch...


                kwargs = rf.meta
                patch_affine = (affine * Affine.translation(cx, cy) *
                                Affine.rotation(angle=-theta) * Affine.translation(-patch_width / 2., -patch_height / 2.))
<commit_msg>Fixed the way I handle 3-pt features -- still no tests<commit_after>""" Extract Patches

This scripts reads in a shapefile and a raster image extracts patches to be used for machine learning algorithms.
If the shapefile consists of points, then a rectangular patch centered at each point is extracted. If the shapefile
consists of lines, then the first and second point are used to determine an orientation, and a patch is extracted
centered at the midpoint of the line and rotated so that the line would be horizontal, pointing in the +x direction.

The ouput folder will always contain a CSV file with the shapefile and raster name, as well as index of the feature in
the shapefile, the latitude and longitude, and image row and column of each feature.


Usage:
    extract_patches.py (-i <shapefile>)...  -r <raster> [--odir <outdir>] [options]
    extract_patches.py -h

Options:
    -i <shapefile>           A shapefile with features indicated. Wildcards are allowed.
    -r <raster>              A  raster source image to extract features from. No wildcards -- if you have multiple
                             raster files consider using `gdalbuildvrt` first to create a virtual mosaic.
    --odir <outdir>          The output folder. [default: ./patches/]
    --otif                   Outout TIF files [default]
    --ojpg                   Output JPG files
    --size <width>,<height>  The size of each extracted patch, in pixels. Please make sure there are no space characters
                             between the arguments unless you put it in quotes! [default: 64,64]
    --scale <scale>          The amount to scale each image by before extracting the patches (e.g.  if the raster is
                             in a higher resolution you may choose to scale it down). [default: 1.0]
    --csv <csvfile>          Create a CSV file with the center lat/lon, center pixel xy, angle, and the patch filename
                             as well as the name of the shapefile source and the raster spource image.
    --silent                 Suppress progress updates.
    --vclip <min,max>        Clip input values to `min` and `max`. Mainly useful when floating point output is not an
                             option (e.g. JPEG output).
    --vstretch <min,max>     Stretch output values to 'min' and 'max', after clipping. Mainly useful when floating point
                             output is not possible. Note that JPEG output must be in 0 to 1.

"""
import hashlib
import os

import numpy
import osgeo
import rasterio
from affine import Affine
from math import degrees, atan2, hypot, floor, ceil
from osgeo import ogr, gdalnumeric
from rasterio._io import RasterReader, RasterUpdater
#from skimage.transform import rescale
from scipy.ndimage import rotate, zoom
from skimage.util.dtype import img_as_ubyte

BREAKING_RELEASE = 0  # Big, rare, when you actually remove deprecated code.
NONBREAKING_RELEASE = 0  # Upload new binaries when you increment this, much less frequently then minor
MINOR_UPDATE = 1  # Increment for pullers
EXPERIMENTAL = 1

VERSION = '{}.{}.{}'.format(BREAKING_RELEASE, NONBREAKING_RELEASE, MINOR_UPDATE)
if EXPERIMENTAL:
    VERSION += 'a{}'.format(EXPERIMENTAL)
__version__ = VERSION

from glob import glob
from docopt import docopt
import logging
from progressbar import Percentage, Bar, ETA, ProgressBar


def collect_filenames(wildcard_list):
    logger = logging.getLogger(__name__)
    results = []
    for wildcard in wildcard_list:
        matching_files = glob(wildcard)
        if len(matching_files) is 0:
            logger.warning("No files matching input specification '{}'".format(wildcard))
        for filename in matching_files:
            results.append(filename)
    return results


def main():
    args = docopt(__doc__, version=VERSION)
    logging.basicConfig(level=logging.INFO)
    logger = logging.getLogger('extract_patches')
    logger.debug('input \n {}'.format(args))

    assert isinstance(logger, logging.Logger)

    shapefiles = collect_filenames(args['-i'])
    if len(shapefiles) == 0:
        logger.error('No matching shapefiles for inoput `{}`'.format(args['-i']))
        return

    raster = args['-r']

    try:
        size = [int(x) for x in args['--size'].split(',')]
        patch_width, patch_height = size
    except:
        logger.error("Unable to parse option '--size'")
        return

    try:
        scale = float(args['--scale'])
        assert scale > 0
    except:
        logger.error("Unable to parse option '--scale'")
        return

    silent = args['--silent']

    output_folder = args['--odir']
    try:
        if not os.path.isdir(output_folder):
            os.makedirs(output_folder)
    except:
        logger.error("Unable to find or create output directory `{}`".format(output_folder))
        return

    if args['--ojpg']:
        fmt = '.jpg'
    else:  # args['--otif']  (default)
        fmt = '.tif'

    clip  = args['--vclip'] is not None
    if clip:
        clipmin, clipmax = [float(x) for x in args['--vclip'].split(',')]

    stretch = args['--vstretch'] is not None
    if stretch:
        stretchmin, stretchmax = [float(x) for x in args['--vstretch'].split(',')]


    # Estimate number of shape features
    count = 0
    if not silent:
        pbar = ProgressBar(len(shapefiles), ['Counting Features:', Percentage(), ' ', Bar(), ' ', ETA()])
        pbar.start()
    for i, s in enumerate(shapefiles):
        vector = ogr.Open(s)
        layer = vector.GetLayer()
        count += layer.GetFeatureCount()
        if not silent:
            pbar.update(i)
    if not silent:
        pbar.finish()

    with rasterio.open(raster) as rf:
        assert isinstance(rf, RasterReader)
        affine = rf.affine
        geo_to_pixels = ~affine

        if not silent:
            pbar = ProgressBar(count, ['Exporting Patches:', Percentage(), ' ', Bar(), ' ', ETA()])
            pbar.start()
        for sf in shapefiles:
            vector = ogr.Open(sf)
            layer = vector.GetLayer()
            for f in layer:
                if not silent:
                    pbar.update(pbar.currval+1)
                geom = f.GetGeometryRef()
                points = geom.GetPoints()
                source = points[0]
                target = points[-1]
                sx, sy = geo_to_pixels * source
                tx, ty = geo_to_pixels * target
                if len(points) == 2:
                    cx, cy = (sx + tx) / 2, (sy + ty) / 2
                else:
                    cx, cy = geo_to_pixels * points[1]
                dx, dy = (tx - sx), (ty - sy)
                theta = degrees(atan2(dy, dx))  # In PIXELS, CCW from +x. Not necessarily CCW from E (or CW from N)
                r1 = hypot(tx - cx, ty - cy)
                r2 = hypot(cx - sx, cy - sy)
                r1, r2 = max(r1, r2), min(r1, r2)  # For 3 points, we assume two radii. Else these are duplicates.
                gx, gy = affine * (cx, cy)  # Geographic coordinates (e.g. lat lon) of the center.

                # We read a square slightly larger than the scaled version of our patch, so that
                # we can safely rotate the raster without missing pixels in the corners.

                box_radius = hypot(patch_width, patch_height) / (2.0 * scale)
                x0, x1 = int(floor(cx - box_radius)), int(ceil(cx + box_radius))
                y0, y1 = int(floor(cy - box_radius)), int(ceil(cy + box_radius))

                                # save patch...


                kwargs = rf.meta
                patch_affine = (affine * Affine.translation(cx, cy) *
                                Affine.rotation(angle=-theta) * Affine.translation(-patch_width / 2., -patch_height / 2.))

{'fraction_non_alphanumeric': 0.09707351891506068, 'fraction_numerical': 0.005515540847446629, 'mean_word_length': 2.1569029086439984, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 29, '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}}
{'dir': 'git-commits-cleaned', 'id': '5479303', 'n_tokens_mistral': 4636, 'n_tokens_neox': 4321, 'n_words': 1616}