expand_point_region.py 9.15 KB
#!/usr/bin/python
# -*- coding: utf-8 -*-

import sys 
import os
import os.path as op
import subprocess
import ogr
import numpy as np
# -----------------------------------------------------------------
# --------- For a polygon perspective -----------------------------
# -----------------------------------------------------------------
def create_buffer_zone(in_shp, out_shp, buffer_dist):
    ''' Create a buffer zone around all the valid points
    A zone is a polygon, which radius is defined by buffer_dist
    
    NOT USED ANYMORE
    '''
    # Get a Layer's Extent
    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.wkbPolygon)

    # Add a class field
    idField = ogr.FieldDefn("class", ogr.OFTInteger)
    outLayer.CreateField(idField)

    # Create the feature and set values
    for point in inLayer:
        ingeom = point.GetGeometryRef()
        geomBuffer = ingeom.Buffer(buffer_dist)
        current_class = point.GetField("class")

        # if class is None, there was an error at the begining,
        # so we discard the point
        if current_class != None:
            featureDefn = outLayer.GetLayerDefn()
            feature = ogr.Feature(featureDefn)
            feature.SetGeometry(geomBuffer)
            feature.SetField("class", current_class)
            outLayer.CreateFeature(feature)

    # Close DataSource
    inDataSource.Destroy()
    outDataSource.Destroy()    
    
    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")
    
# -----------------------------------------------------------------
# --------- For a points perspective ------------------------------
# -----------------------------------------------------------------

def km_to_deg_convert(km, direction, latitude=0):
    ''' Convert roughly a distance in km to degrees
    The latitude parameter should be in degrees, and is needed for the 
    computation in case of longitude
    '''
    if direction == 'lat':
        deg = km/110.574
    elif direction == 'lon':
        deg = km/(111.320*np.cos(np.radians(latitude)))
    else:
        deg = 0
    return deg
    #~ Latitude: 1 deg = 110.574 km
    #~ Longitude: 1 deg = 111.320*cos(latitude) km
    


def create_neighbors(in_shp, out_shp, max_dist_X, max_dist_Y, nb_samples):
    ''' Create a neighbourhoud around all the points in a shapefile
    For each point, a grid around it is created, from -max_dist to +max_dist
    in each direction, with nb_samples in each direction
    
    NOT USED ANYMORE
    '''
    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")
    
    # Watch out : sometimes the unit is meter, sometimes degree
    srs_unit = srs.GetAttrValue('unit')
    
    if srs_unit == 'Meter':
        print('Unit is meter, no change')
    elif srs_unit == 'Degree':
        print('Unit is degree, needs to be converted')

    # 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)

    print('{} neighbors will be created'.format(nb_samples*nb_samples*len(inLayer)) )

    # Create the feature and set values
    for point in inLayer:
        current_class = point.GetField("class")
        # if class is None, there was an error at the begining,
        # so we discard the point
        if current_class != None:
            ingeom = point.GetGeometryRef()
            
            # define the extent of the grid
            X_range = np.linspace(start=ingeom.GetX(0)-max_dist_X, stop=ingeom.GetX(0)+max_dist_X, num=nb_samples)
            Y_range = np.linspace(start=ingeom.GetY(0)-max_dist_Y, stop=ingeom.GetY(0)+max_dist_Y, num=nb_samples)
            
            # create the grid
            for i in range(len(X_range)):
                for j in range(len(X_range)):
                    outpoint = ogr.Geometry(ogr.wkbPoint)
                    outpoint.AddPoint(X_range[i], Y_range[j])
                    
                    featureDefn = outLayer.GetLayerDefn()
                    feature = ogr.Feature(featureDefn)
                    feature.SetGeometry(outpoint)
                    feature.SetField("class", current_class)
                    outLayer.CreateFeature(feature)

    # Close DataSource
    inDataSource.Destroy()
    outDataSource.Destroy()    
    
    return


def create_squares(in_shp, out_shp, max_dist_X, max_dist_Y):
    ''' 
    Create a neighbourhoud around all the points in a shapefile
    For each point, a square around it is created, from -max_dist to +max_dist
    in each direction
    '''
    inDriver = ogr.GetDriverByName("ESRI Shapefile")
    inDataSource = inDriver.Open(in_shp, 0)
    inLayer = inDataSource.GetLayer()
    
    srs = inLayer.GetSpatialRef()
    
    # Watch out : sometimes the unit is meter, sometimes degree
    srs_unit = srs.GetAttrValue('unit')
    
    if srs_unit == 'Meter':
        print('Unit is meter, no change')
    elif srs_unit == 'Degree':
        print('Unit is degree, needs to be converted')

    # 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.wkbPolygon)

    # Add a class field
    classField = ogr.FieldDefn("class", ogr.OFTInteger)
    outLayer.CreateField(classField)
    print('{} squares will be created'.format(len(inLayer)) )
    
    
    # Create the feature and set values
    for point in inLayer:
        current_class = point.GetField("class")
        if current_class != None:
            ingeom = point.GetGeometryRef()
            
            Xpoint = ingeom.GetX(0)
            Ypoint = ingeom.GetY(0)
            
            # definition of the square sides
            left = Xpoint-max_dist_X
            right = Xpoint + max_dist_X
            top = Ypoint+max_dist_Y
            bottom = Ypoint-max_dist_Y
            
            border = ogr.Geometry(ogr.wkbLinearRing)
            # 4 corners, needs to be closed with the 5th point
            border.AddPoint(left, top)
            border.AddPoint(right, top)
            border.AddPoint(right, bottom)
            border.AddPoint(left, bottom)
            border.AddPoint(left, top)
            
            # create the polygon
            poly = ogr.Geometry(ogr.wkbPolygon)
            poly.AddGeometry(border)
            
            # assign the class
            featureDefn = outLayer.GetLayerDefn()
            feature = ogr.Feature(featureDefn)
            feature.SetGeometry(poly)
            feature.SetField("class", current_class)
            outLayer.CreateFeature(feature)    
        

    # Close DataSource
    inDataSource.Destroy()
    outDataSource.Destroy()    
    
    return





def main():
    in_shp = '/mnt/data/home/baetensl/OTB_codes/OTB_commands/Full_orleans/In_data/Masks/human_reference.shp'
    out_shp = '/mnt/data/home/baetensl/OTB_codes/OTB_commands/Full_orleans/In_data/Masks/NEW_delete.shp'
    
    in_shp = '/mnt/data/home/baetensl/classification_clouds/Data/Orleans_all/In_data/Masks/land.shp'
    out_shp = '/mnt/data/home/baetensl/classification_clouds/Data/Orleans_all/In_data/Masks/SQUARES.shp'
    

    #~ in_shp = '/mnt/data/home/baetensl/classification_clouds/Data/Orleans_all/Intermediate/train_points.shp'
    #~ out_shp = '/mnt/data/home/baetensl/classification_clouds/Data/Orleans_all/Intermediate/train_points222.shp'

    #~ createBuffer(in_shp, out_shp, buffer_dist = 0.01)
    #~ translate_points(in_shp, out_shp, translation_dist = 0.01)
    
    km = 50000
    degX = km_to_deg_convert(km, direction='lon', latitude=48)
    degY = km_to_deg_convert(km, direction='lat', latitude=48)

    degX = 100
    degY = 100
    #~ create_neighbors(in_shp, out_shp, max_dist_X = degX, max_dist_Y = degY, nb_samples=5)
    create_squares(in_shp, out_shp, max_dist_X = degX, max_dist_Y = degY)
    
if __name__ == '__main__':
    main()