3

I am trying to create some polygons where some of the sections correspond to curves/arcs. I have been provided with start point, end point, centre point, radius and direction (clockwise/ anti-clockwise) to create these curved sections. See below for reference:

enter image description here

I have seen this post that has helped massively but unfortunately I do not know how to derive my arc midpoint from the values I have been given (centre point and radius). I would like to do this programmatically using Python/PyQGIS but somehow, I am unable to find a solution. So far, I have managed to manually trace the midpoint section for the curve and used it to create my lineString layer.

Is there a pythonic way to achieve this and convert to polygon?

This is my code so far:

# Import modules
from qgis.core import *
from qgis.utils import *

# coord for curved section
start_lat, start_lon = 39.08626748, -2.66563606
end_lat, end_lon = 39.20366816, -2.51489454
# This mid point has been traced manually to recreate the curve
mid_curv_lat, mid_curv_lon = 39.20366816, -2.65091330

# coord for squared section
point3_lat, point3_lon = 39.11293834, -2.41907255
point4_lat, point4_lon = 38.98554723, -2.55822285

# coord and radius
centre_lat, centre_lon = 39.13921171, -2.58290392
radius_km = 9.26

    # Begin a new project
    iface.newProject()
    
    # Create a QgsCircularString
    circularRing = QgsCircularString()
    # Set first point, intermediate point for curvature and end point
    circularRing.setPoints([
        QgsPoint(start_lon, start_lat),
        QgsPoint(mid_curv_lon, mid_curv_lat),
        QgsPoint(end_lon, end_lat)]
    )
    
    # Create geometry using the instance of QgsCircularString
    curve = QgsGeometry(circularRing)
    
    # Create a QgsLineString
    lineString = QgsGeometry.fromPolyline([QgsPoint(end_lon, end_lat), 
                 QgsPoint(point3_lon, point3_lat),
                 QgsPoint(point4_lon, point4_lat),
                 QgsPoint(start_lon, start_lat)]
    )
    
    # Create curve and line features
    fet_curve = QgsFeature()
    fet_curve.setGeometry(curve)
    
    fet_line = QgsFeature()
    fet_line.setGeometry(lineString)
    
    # Create a memory layer
    layer = QgsVectorLayer(
        "LineString?crs=epsg:4326&field=id:integer&field=name:string(20)&index=yes",
        "temporary_points",
        "memory"
    )
    # Add the layer
    QgsProject.instance().addMapLayer(layer)
    # Add the feature to the layer provider
    pr = layer.dataProvider()
    pr.addFeatures([fet_curve, fet_line])
    
    # Update extent
    layer.updateExtents()
    
    # Zoom to extent
    iface.mapCanvas().setExtent(layer.extent())
    iface.mapCanvas().refresh()

EDIT: data sample

|Label|Position|Order|CoordPair         |Orientation|Radius_Km|
|-----|--------|-----|------------------|-----------|---------|
|C    |        |1    |39.08627, -2.66564|           |         |
|C    |        |2    |39.20367, -2.51489|           |         |
|C    |        |3    |39.11294, -2.41907|           |         |
|C    |        |4    |38.98555, -2.55822|           |         |
|C    |Centre  |     |39.13921, -2.5829 |1          |9.26     |
|A    |        |5    |38.91119, -2.25953|           |         |
|A    |        |1    |38.99365, -2.21096|           |         |
|A    |        |2    |38.94734, -2.09687|           |         |
|A    |        |3    |38.89763, -2.10139|           |         |
|A    |        |4    |38.86149, -2.18385|           |         |
|A    |Centre  |     |38.927, -2.11607  |1          |9.26     |
|B    |        |2    |38.74627, -2.41767|           |         |
|B    |        |3    |38.74062, -2.26179|           |         |
|B    |        |4    |38.68866, -2.25953|           |         |
|B    |        |1    |38.69769, -2.42445|           |         |
|B    |Centre1 |     |38.72367, -2.44253|1          |9.26     |
|B    |Centre2 |     |38.71464, -2.25388|1          |9.26     |
6
  • What if you separately draw the circular and rectangular portions, then get the union of those geometries?
    – Kartograaf
    Jan 5 at 19:10
  • 1
    How did you calculate radius (km) from coordinates with degree unit? Jan 5 at 20:03
  • @Kartograaf that's not feasible as some of the "rectangular" polygons are not regular and will create unwanted sections. Jan 6 at 11:54
  • @KadirŞahbaz I manually traced the mid point of the arc section after creating a buffer around the point for illustrative purposes only. I have not found a way to create the arc with the values only. Jan 6 at 11:55
  • Can you share any sample data?
    – Kartograaf
    Jan 6 at 18:52
2

I believe QgsCircularString.fromTwoPointsAndCenter() is a method that could help you construct those arcs. It creates a circular string consisting of a single arc which represents the curve from one point to another with a specified center point.

I've updated your example using the mentioned method and also addressed your problem of converting the exterior line to a polygon.

# coord for curved section
start_lat, start_lon = 39.08626748, -2.66563606
end_lat, end_lon = 39.20366816, -2.51489454

# coord for squared section
point3_lat, point3_lon = 39.11293834, -2.41907255
point4_lat, point4_lon = 38.98554723, -2.55822285

# center coord
centre_lat, centre_lon = 39.13921171, -2.58290392

# Create the arc segment
arc = QgsCircularString.fromTwoPointsAndCenter(
    QgsPoint(start_lon, start_lat),
    QgsPoint(end_lon, end_lat),
    QgsPoint(centre_lon, centre_lat)
)

# Creates the polygon segment
line = QgsLineString(
    [QgsPoint(end_lon, end_lat), QgsPoint(point3_lon, point3_lat),
    QgsPoint(point4_lon, point4_lat), QgsPoint(start_lon, start_lat)]
)

# join both segments together
curve = line.toCurveType()
curve.addCurve(arc)

# create polygon geometry from exterior line
polygon = QgsPolygon()
polygon.setExteriorRing(curve)
geom = QgsGeometry(polygon)

# Create polygon feature
poly_feat = QgsFeature()
poly_feat.setGeometry(polygon)

# Create a memory layer
layer = QgsVectorLayer(
    "Polygon?crs=epsg:4326&field=id:integer&field=name:string(20)&index=yes",
    "temporary_points",
    "memory"
)
# Add the layer
QgsProject.instance().addMapLayer(layer)
# Add the feature to the layer provider
pr = layer.dataProvider()
pr.addFeatures([poly_feat])

# Update extent
layer.updateExtents()

# Zoom to extent
iface.mapCanvas().setExtent(layer.extent())
iface.mapCanvas().refresh()

arcs

3
  • thanks for your response but unfortunately I need 9.26km from the centre at any point in the curve. If you create a buffer of 9.26km from the centre and overlay the result with the output from your script you can see where this is not happening. Is there any other method that perhaps will achieve this? Jan 20 at 11:51
  • 1
    I really can't fully explain why that doesn't work, but it has something to do with the projection and that EPSG:4326 uses degrees as map units instead of meters. So I recommend that you use a local projection which uses meters as unit, e.g. UTM zone 30N (EPSG:25830). Running the script again with projected coordinates does output what you are trying to achieve. Take a look at this image where I marked the buffer: imgur.com/a/O3XYrBG
    – CodeBard
    Jan 20 at 16:52
  • if I transform the coordinates your solution works. Thanks for this. It's much appreciated. Jan 24 at 9:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.