"""
Clip CanopyHeight data to the same extent/resolution/geo-coordinates of the flightline DSM data
The CanopyHeight data are later used to model the offsets of the tgt img relative to the reference img.
"""
import gdal,glob, sys
import os
import time
import numpy as np, pandas as pd
from os.path import expanduser

Dir_in = r'Z:\townsenduser-rw\HyspexPro\Output\Cheesehead_V2'
sub_ls = glob.glob(Dir_in+'/CHEESE*/')
# remove -1_20190711
sub_ls = [x for x in sub_ls if '-1_20190711' not in x]
sub_ls.sort()
CH = r'Z:\townsenduser-rw\projects\CHEESEHEAD\CH_Lidar' + '/CH_CanopyHeight_30km_1m'
#------ Step 1: Convert DSM extent to shape file-------------------------------
for sub in sub_ls[3:]:
    print('process subfolder %s' %sub)
    imgs = glob.glob(sub + 'Merge/*_DEM')
    imgs.sort()
    for img in imgs:
        outname = img.replace('_DEM','_DEM_shp.shp')
        CMD = 'gdaltindex ' + outname + ' ' + img
        os.system(CMD)

        #------- Step 2: Clip the canopy height data based on the shape files
        CH_name = img.replace('_DEM', '_CanopyHeight.tif')
        # Get the DEM resolution
        DEM = gdal.Open(img)
        gt = DEM.GetGeoTransform()
        xres = abs(gt[1])
        yres = abs(gt[5])
        del DEM
        clip_CMD = 'gdalwarp ' + '-tr ' + str(xres) + ' ' + str(yres) + ' -cutline ' + outname + ' -crop_to_cutline ' + CH + ' ' + CH_name
        os.system(clip_CMD)