perp.py

def remove_layers(layer=None):
    qinst=QgsProject.instance()

    lyrs=['perp', 'fwlshift', 'slope_vector'] if not layer else [layer]
    for name in lyrs:
        try:
            for l in qinst.mapLayersByName(name):
                qinst.removeMapLayer(l.id())
        except:
            print("Layer not found (not an error)")
            


def get_centroid(island_area_poly='island_centroid'):
    qinst=QgsProject.instance()
    # make sure all layers are projected to the projects CRS (If not just export them with projects CRS and reload)
    l1=qinst.mapLayersByName(island_area_poly)[0]
    l1.selectByExpression('"ID"=1 OR "FID"=1')
    sel1=l1.selectedFeatures()
    if(len(sel1)!=1):
        print("ERROR we need exactly 1 feature with ID=1 (data {}".format(str(sel1)))
        raise Exception
    centroid_point = sel1[0].geometry().asPoint()
    return centroid_point


def get_100mline(bathy_contour_parallel_line='100m_temp'):
    remove_layers(layer=bathy_contour_parallel_line)
    qinst=QgsProject.instance()
    l2=qinst.mapLayersByName(bathy_contour_parallel_line)[0]
    l2.selectByExpression('"ID"=1 OR "FID"=1')
    sel2=l2.selectedFeatures()
    if(len(sel2)!=1):
        print("ERROR we need exactly 1 feature with ID=1 (data {}".format(str(sel2)))
        raise Exception
    feat_line = sel2[0]
    return feat_line

def get_2mcentroid(fwlcont_centroid='fwlcentroid'):
    qinst=QgsProject.instance()
    l3=qinst.mapLayersByName(fwlcont_centroid)[0]
    l3.selectByExpression('"ID"=1 OR "FID"=1')
    sel3=l3.selectedFeatures()
    if(len(sel3)!=1):
        print("ERROR we need exactly 1 feature with ID=1 (data {}".format(str(sel3)))
        raise Exception
    fwl_point=sel3[0]
    return fwl_point.geometry().asPoint()

def get_gradient_vector_layer(vector_layer_name='bathymetry_maldives_gradient_vectors2'):
    qinst=QgsProject.instance()
    try:
        return qinst.mapLayersByName(vector_layer_name)[0]
    except:
        print (vector_layer_name, "not found in the project")
        raise Exception



def shift_based_on_100m_line_perp(feat_line, centroid_point, style="../arrow_style_for_100m_bathy.qml"):
    qinst=QgsProject.instance()
    geom = feat_line.geometry().closestSegmentWithContext(centroid_point.geometry().asPoint())
    #geom is a tupla. I need only second term, geom[1], that is a QgsPoint
    closest_line = QgsGeometry.fromPolyline([ QgsPoint(centroid_point.geometry().asPoint()), QgsPoint(geom[1]) ])
    #print(closest_line)
    crs=qinst.crs().authid() # get project's crs
    capa = QgsVectorLayer("LineString?crs={}".format(crs), "perp", "memory")

    pr = capa.dataProvider()
    pr.addAttributes([
                      QgsField("id",  QVariant.Int),
                      QgsField("name", QVariant.String)])
    capa.updateFields()

    # Add the layer in QGIS project
    qinst.addMapLayer(capa)

    # Start editing layer
    capa.startEditing()
    feat = QgsFeature(capa.fields()) # Create the feature
    feat.setAttribute("id", 1) # set attributes
    feat.setAttribute("name", "perp")
    feat.setGeometry(closest_line)
    capa.addFeature(feat) # add the feature to the layer
    capa.endEditCommand() # Stop editing
    capa.commitChanges() # Save changes

    capa.loadNamedStyle(qinst.readPath("/")+"/"+style)
    return qinst



def draw_vector_sum_direction(centroid_point, vector_layer, layer="slope_vector", style="../arrow_style_slope_vector2.qml", radius=10):
    qinst=QgsProject.instance()    
    remove_layers(layer)
    # 
    RAD=radius*1000
    geom_buffer = QgsGeometry.fromPointXY(centroid_point).buffer(RAD, -1)
    # remove past selections
    vector_layer.removeSelection()
    for feature in vector_layer.getFeatures():
        if feature.geometry().intersects(geom_buffer):
            #print(feature.id())
            vector_layer.select(feature.id())
    vector_layer.updateFields()


    # loop again - inefficent could have been done in the previous loop, but doing to keep modularity
    SEX=0
    SEY=0
    NI =len(vector_layer.selectedFeatures())
    for feature in vector_layer.selectedFeatures():
        SEX+=feature['EX']
        SEY+=feature['EY']

    
    crs=qinst.crs().authid() # get project's crs
    v_layer = QgsVectorLayer("LineString?crs={}".format(crs), layer, "memory")
    pr = v_layer.dataProvider()
    seg = QgsFeature()
    endpoint=centroid_point+QgsVector(SEX*10,SEY*10)
    seg.setGeometry(
        QgsGeometry.fromPolyline([QgsPoint(centroid_point), QgsPoint(endpoint)] ))
    pr.addFeatures([ seg ])
    qinst.addMapLayers([v_layer])
    v_layer.loadNamedStyle(qinst.readPath("/")+"/"+style)
    return endpoint


def draw_centroid_shift(centroid_point, fwl_point, layer="fwlshift", style="../arrow_style_for_centroid_shift.qml"):
    qinst=QgsProject.instance()   
    remove_layers(layer)
    crs=qinst.crs().authid() # get project's crs
    # draw the centroild shift

    v_layer = QgsVectorLayer("LineString?crs={}".format(crs),  layer, "memory")
    pr = v_layer.dataProvider()
    seg = QgsFeature()
    seg.setGeometry(
        QgsGeometry.fromPolyline([QgsPoint(centroid_point), 
                                  QgsPoint(fwl_point)
                                  ]
                                 )
        )
    pr.addFeatures([ seg ])
    qinst.addMapLayers([v_layer])
    v_layer.loadNamedStyle(qinst.readPath("/")+"/"+style)    

def mark_angle(centroid_point, fwl_point, vector_point):
    azimuth1=centroid_point.azimuth(fwl_point)
    azimuth2=centroid_point.azimuth(vector_point)
    diff=abs(azimuth1-azimuth2)
    if diff > 180:
        diff=360.-diff
    return diff
    

#remove_layers()
centroid_point = get_centroid(island_area_poly='island_centroid')
vector_layer=get_gradient_vector_layer(vector_layer_name='gradient_vectors')
fwl_point = get_2mcentroid(fwlcont_centroid='fwlcentroid')
#perp_line = get_100mline(bathy_contour_parallel_line='100m_temp')
#shift_based_on_100m_line_perp(perp_line, centroid_point, style="../arrow_style_for_100m_bathy.qml")
draw_centroid_shift(centroid_point, fwl_point, layer="fwlshift", style="../arrow_style_for_centroid_shift.qml")
radius=10 # km
vector_point=draw_vector_sum_direction(centroid_point, vector_layer, radius=radius, layer="slope_vector", style="../arrow_style_slope_vector2.qml", )
mark_angle(centroid_point, fwl_point, vector_point)