# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python [conda env:dsi]
#     language: python
#     name: conda-env-dsi-py
# ---

# +
import json
from os import path, walk, makedirs
from sys import exit, stderr

from cv2 import fillPoly, imwrite
import numpy as np
from shapely import wkt
from shapely.geometry import mapping, Polygon
from skimage.io import imread
from tqdm import tqdm
import imantics 

# This removes the massive amount of scikit warnings of "low contrast images"
import warnings
warnings.filterwarnings("ignore", category=UserWarning)


# -

# ### intersection with images: crop and generate masks

# This section crops the geotiffs according to the bounding box for each area of interest. It then finds the intersection of the damaged building polygons with the satellite image and creates a pixel mask of 0s and 1s, with 1s indicating the pixel is part of a damaged building.

# + jupyter={"outputs_hidden": false}
def crop_images(coordinates, image_dir_path, img_out_path, series):

    img_list = [os.path.join(image_dir_path, image_name) 
                   for image_name in os.listdir(image_dir_path) 
                   if image_name.lower().endswith('tif')]

    image_names = [image_name[:-4] for image_name in os.listdir(image_dir_path) if image_name.lower().endswith('tif')]

    i = 0
    for image in img_list:
        geom = [{'type': 'Polygon', 'coordinates': [coordinates]}]

        with rasterio.open(image) as src:
            out_image, out_transform = mask(src, geom, crop = True)
        out_meta = src.meta.copy()

        # save the resulting raster  
        out_meta.update({"driver": "GTiff",
                        "height": out_image.shape[1],
                        "width": out_image.shape[2],
                        "transform": out_transform})

        with rasterio.open(img_out_path + image_names[i] + '_crop_' + str(series) + '.tif', "w", **out_meta) as dest:
            dest.write(out_image)  
        i +=1


# -

def polycoords(poly):
    """Convert a polygon into the format expected by OpenCV
    """
    if poly.type in ['MultiPolygon', 'GeometryCollection']:
        return [np.array(p.exterior.coords) for p in poly if p.type == 'Polygon']
    elif poly.type == 'Polygon':
        return [np.array(poly.exterior.coords)]
    else:
        print('Encountered unrecognized geometry type {}. Ignoring.'.format(poly.type))
        return []


def generate_masks(image_path, building_polygon_df):
    
    def make_mask_from_many(img_shape, poly_list):
        """Make a mask from multiple polygons"""
        poly_pts = [polycoords(poly) for poly in coord_list if poly]
        # note: deal with rounding here??
        polys = [[x.astype(int) for x in j] for j in poly_pts]
        polys = [np.asarray(poly) for poly in polys]
        # Create an empty mask and then fill in the polygons
        mask = np.zeros(img_shape[:2])
        cv2.fillPoly(mask, polys, 1)
        return mask.astype('uint8')

    # create mask for given image

    src = rasterio.open(image_path)
    img = src.read().transpose([1,2,0])
    img_bounds = shapely.geometry.box(*src.bounds)
    img_transform = list(np.array(~src.affine)[[0, 1, 3, 4, 2, 5]])

    coord_list = []
    
    for poly in collapsed.geometry:

        mask_poly = poly.intersection(img_bounds)
        mask_poly_pxcoords = shapely.affinity.affine_transform(mask_poly, img_transform)
        coord_list.append(mask_poly_pxcoords)
            
    mask = make_mask_from_many(img.shape, coord_list)
    # removing first row which is a border (remnant of cropping)
    mask = mask[1:, :]

    return mask


# + jupyter={"outputs_hidden": false}
# crop talcoban city
coords = list(talcoban.boundary.coords)

images_in = '/mnt/merged_images/'
imgages_out = '/mnt/cropped_images/'
crop_images(coords, images_in, imgages_out, 1)

# + jupyter={"outputs_hidden": false}
# print(talcoban.centroid)

# + jupyter={"outputs_hidden": false}
# samar
coords = list(samar.boundary.coords)

images_in = '/mnt/merged_images/'
imgages_out = '/mnt/cropped_images/'
crop_images(coords, images_in, imgages_out, 2)

# +
# palo
coords = list(palo.boundary.coords)

images_in = '/mnt/merged_images/'
imgages_out = '/mnt/cropped_images/'
crop_images(coords, images_in, imgages_out, 3)

# + jupyter={"outputs_hidden": false}
print(palo.centroid)

# +
# guiuan
coords = list(guiuan.boundary.coords)

images_in = '/mnt/merged_images/'
imgages_out = '/mnt/cropped_images/'
crop_images(coords, images_in, imgages_out, 4)

# + jupyter={"outputs_hidden": false}
# generate masks from cropped sharpened post image
image_in = '/mnt/cropped_images/POST_pan_sharp_crop_1.tif'
mask1 = generate_masks(image_in, collapsed)
image_in = '/mnt/cropped_images/POST_pan_sharp_crop_2.tif'
mask2 = generate_masks(image_in, collapsed)
image_in = '/mnt/cropped_images/POST_pan_sharp_crop_3.tif'
mask3 = generate_masks(image_in, collapsed)
image_in = '/mnt/cropped_images/POST_pan_sharp_crop_4.tif'
mask4 = generate_masks(image_in, collapsed)

# + jupyter={"outputs_hidden": true}

# -

# **a few checks:**

# Checks confirm mask shapes are reasonable and that masks have nonzero values.

# + jupyter={"outputs_hidden": false}
mask1.shape

# + jupyter={"outputs_hidden": false}
np.nonzero(mask1)

# + jupyter={"outputs_hidden": false}
mask2.shape

# + jupyter={"outputs_hidden": false}
# np.nonzero(mask2)

# + jupyter={"outputs_hidden": false}
mask3.shape

# + jupyter={"outputs_hidden": false}
np.nonzero(mask3)

# + jupyter={"outputs_hidden": false}
mask4.shape

# + jupyter={"outputs_hidden": false}
np.nonzero(mask4)

# + jupyter={"outputs_hidden": false}
# percent of nonzeros for mask1
len(np.nonzero(mask1)[0]) / (len(mask1[0])*len(mask1[1])) * 100

# + jupyter={"outputs_hidden": false}
# percent of nonzeros for mask2
len(np.nonzero(mask2)[0]) / (len(mask2[0])*len(mask2[1])) * 100

# + jupyter={"outputs_hidden": false}
# percent of nonzeros for mask3
len(np.nonzero(mask3)[0]) / (len(mask3[0])*len(mask3[1])) * 100

# + jupyter={"outputs_hidden": false}
# percent of nonzeros for mask4
len(np.nonzero(mask4)[0]) / (len(mask4[0])*len(mask4[1])) * 100

StarcoderdataJupyterScriptsDedupFiltered