layers_creation.py 7.47 KB
#!/usr/bin/python
# -*- coding: utf-8 -*-

 
import os
import os.path as op
import ogr
from osgeo import gdal,osr
import json
import otbApplication
import L1C_band_composition
import subprocess
import tempfile

def empty_shapefile_creation(in_tif, out_shp_list, geometry_type = 'point'):
    ''' 
    Create empty shapefiles based on the SRS of the in_tif
    '''

    # get the projection and SRS
    ds=gdal.Open(r'{}'.format(in_tif))
    prj=ds.GetProjection()
    srs=osr.SpatialReference(wkt=prj)

    for k in range(len(out_shp_list)):
        out_shp = out_shp_list[k]

        outDriver = ogr.GetDriverByName("ESRI Shapefile")
        # Remove output shapefile if it already exists
        if os.path.exists(out_shp):
            outDriver.DeleteDataSource(out_shp)

        # Create the output shapefile
        outDataSource = outDriver.CreateDataSource(out_shp)
        if geometry_type == 'point':
            outLayer = outDataSource.CreateLayer("buff_layer", srs, geom_type=ogr.wkbPoint)
        elif geometry_type == 'polygon':
            outLayer = outDataSource.CreateLayer("buff_layer", srs, geom_type=ogr.wkbPolygon)
        else:
            print('Please enter a valid geometry type')
        # Add a class field
        classField = ogr.FieldDefn("class", ogr.OFTInteger)
        outLayer.CreateField(classField)
        
        # Close DataSource
        outDataSource.Destroy()    
    
    return
    
def create_all_classes_empty_layers(global_parameters, force = False):
    '''
    Create all the empty layers files
    '''
    print('  Creation of the empty layers')
    main_dir = global_parameters["user_choices"]["main_dir"]
    
    # append from the configuration file the masks to create
    layers_to_create = []
    for mask_name, mask_values in global_parameters["masks"].iteritems():
        layers_to_create.append(op.join(main_dir, 'In_data', 'Masks', mask_values["shp"]))
        
    in_tif = op.join(main_dir, 'In_data', 'Image', global_parameters["user_choices"]["raw_img"])
    
    # check if they already exist
    if any([op.exists(f) for f in layers_to_create]) and force == False:
        print('Layers already present, use -force to erase and replace')
        return
        
    empty_shapefile_creation(in_tif, layers_to_create)
    
    # same for the no_data, but separate as it is a polygon type
    no_data_layer = [op.join(main_dir, 'In_data', 'Masks', global_parameters["general"]["no_data_mask"])]
    empty_shapefile_creation(in_tif, no_data_layer, geometry_type = 'polygon')
    print('Done')
    
 
def populate_layer(in_shp, out_shp, x_coords, y_coords):
    '''
    Populate a layer with points from their geo coordinates
    '''
    inDriver = ogr.GetDriverByName("ESRI Shapefile")
    inDataSource = inDriver.Open(in_shp, 0)
    inLayer = inDataSource.GetLayer()
    
    srs = inLayer.GetSpatialRef()

    # Save extent to a new Shapefile
    outDriver = ogr.GetDriverByName("ESRI Shapefile")
    # Remove output shapefile if it already exists
    if os.path.exists(out_shp):
        outDriver.DeleteDataSource(out_shp)

    # Create the output shapefile
    outDataSource = outDriver.CreateDataSource(out_shp)
    outLayer = outDataSource.CreateLayer("buff_layer", srs, geom_type=ogr.wkbPoint)

    # Add a class field
    classField = ogr.FieldDefn("class", ogr.OFTInteger)
    outLayer.CreateField(classField)
    
    # Create the feature and set values
    for point_nb in range(len(x_coords)):
        # set the location of the point
        point = ogr.Geometry(ogr.wkbPoint)
        point.SetPoint(0, x_coords[point_nb], y_coords[point_nb])
        
        featureDefn = outLayer.GetLayerDefn()
        feature = ogr.Feature(featureDefn)
        feature.SetGeometry(point)
        outLayer.CreateFeature(feature)

    # Close DataSource
    inDataSource.Destroy()
    outDataSource.Destroy()    
    
    return
    
    
    
def fill_precomputed_classes():
    '''
    DEPRECATED
    '''
    in_tif = '/mnt/data/home/baetensl/clouds_detection_git/Data_ALCD/RailroadValley_11SPC_20170501/In_data/Image/RailroadValley_bands.tif'
    out_shp_list = ['/mnt/data/home/baetensl/clouds_detection_git/Various_data/test_out.shp']
    empty_shapefile_creation(in_tif, out_shp_list, geometry_type = 'point')
    
    in_shp = '/mnt/data/home/baetensl/clouds_detection_git/Various_data/test_out.shp'
    out_shp = '/mnt/data/home/baetensl/clouds_detection_git/Various_data/filled.shp'
    x_coords = [620000, 700000]
    y_coords = [4250000, 4280000]
    populate_layer(in_shp, out_shp, x_coords, y_coords)
                
    
    
    
    
    

    
    
def create_no_data_shp(global_parameters, force = False):    
    '''
    Create automatically polygons overs the no_data pixels
    in both the clear and cloudy date
    '''
    main_dir = global_parameters["user_choices"]["main_dir"]
    
    if op.exists(op.join(main_dir, 'In_data', 'Masks', 'no_data.shp')) and force == False:
        print('No_data layer already present, use -force to erase and replace')
        return

    tmp_name = next(tempfile._get_candidate_names())
    tmp_tif = op.join('tmp', 'no_data_mask_{}.tif'.format(tmp_name))

    # Create the temporary no data TIF
    L1C_band_composition.create_no_data_tif(global_parameters, tmp_tif)
    
    # polygonize the raster
    print("  Polygonization")
    tmp_shp = op.join('tmp', 'no_data_mask_{}.shp'.format(tmp_name))
    out_layer = 'no_data_shape'
    command = 'gdal_polygonize.py {} -mask {} -f "ESRI Shapefile" {} {} class'.format(tmp_tif, tmp_tif, tmp_shp, out_layer)
    subprocess.call(command, shell=True)
    
    # simplify the polygons
    out_shp = op.join(main_dir, 'In_data', 'Masks', 'no_data.shp')
    simplify_geometry(tmp_shp, out_shp, tolerance = 100)
    return

def simplify_geometry(in_shp, out_shp, tolerance = 100):
    ''' Simplification of a shapefile
    This allows to have lighter polygons, enhancing the rapidty
    '''
    print("  Geometry simplification")    
    command = "ogr2ogr {} {} -simplify {}".format(out_shp, in_shp, str(tolerance))
    subprocess.call(command, shell=True)
    print("Done")
    
    
    
    
    
def main():
    global_parameters = json.load(open(op.join('parameters_files','global_parameters.json')))    

    create_all_classes_empty_layers(global_parameters, force = False)
    
    return
    #~ in_tif = '/mnt/data/home/baetensl/classification_clouds/Data/Full_orleans/In_data/Image/orleans_bands.tif'
    create_no_data_shp(global_parameters)
    
    return
    main_dir = '/mnt/data/home/baetensl/clouds_detection_git/Data_ALCD/Arles_31TFJ_20171002'
    in_tif = op.join(main_dir, 'In_data', 'Image', 'Arles_bands.tif')
    out_shp = op.join(main_dir, 'Intermediate', 'no_data_autom.shp')
    
    no_data_preprocess(in_tif, out_shp)
    
    #~ fill_precomputed_classes()
    
    return
    in_tif = '/mnt/data/home/baetensl/clouds_detection_git/Data_ALCD/RailroadValley_11SPC_20170501/In_data/Image/RailroadValley_bands.tif'
    
    layers_to_create = []
    global_parameters = json.load(open(op.join('parameters_files','global_parameters.json')))
    
    create_all_classes_empty_layers(global_parameters)
        
    #~ for key, value in global_parameters["masks"].iteritems():
        #~ layers_to_create.append(op.join(global_parameters["general"]["main_dir"], 'In_data', 'Masks', value))
    #~ print(layers_to_create[0])
    
    #~ out_shp_list = ['/mnt/data/home/baetensl/classification_clouds/Data/Full_orleans/In_data/Image/land_shp.shp']
    #~ empty_shapefile_creation(in_tif, layers_to_create)
    
if __name__ == '__main__':
    main()