0

If two polylines intersect with each other and there are more than 2 intersect points, there is more than 1 closed area will be created.

How can I get these areas as polygons in ArcGIS Desktop or Python? For example, the yellow areas are created after two polylines intersect each other. I need to calculate the area of the yellow polygon.

Area of differences

1
  • I think this question would benefit from the inclusion of a picture that illustrates your input and desired output.
    – PolyGeo
    Aug 1, 2019 at 5:48

2 Answers 2

1

The current solution:

  1. Merge modeled stream network line features into one single feature Shapefile A using the Merge tool;
  2. Merge NHD flowline features (actual flowlines) into one single feature Shapefile B using the Merge tool;
  3. Merge the Shapefile A with the Shapefile B to obtain Shapefile C;
  4. Convert polyline Shapefile C into polygon Shapefile D using the FeatureToPolygon tool; and
  5. Add a new field and calculate the area of the polygon features using the CalculateGeometry tool. The total area of differences is the sum of all polygon features

So far, we haven't found a way to do this in Python GDAL.

2
  • This is overkill simply use the feature to polygon tool.
    – Hornbydd
    Aug 1, 2019 at 15:24
  • Indeed, we simplified the solution using the polyline to polygon tool:
    – Chang
    Jan 15, 2022 at 21:44
1

Polygonize the lines using Shapely

Quoting the following snippets from Shapely manual

shapely.ops.polygonize(lines)
#Returns an iterator over polygons constructed from the input lines.

There is also another function called polygonize_full.

shapely.ops.polygonize_full(lines)

It creates polygons from a source of lines, returning the polygons and leftover geometries.The source may be a MultiLineString, a sequence of LineString objects, or a sequence of objects than can be adapted to LineStrings.Returns a tuple of objects: (polygons, dangles, cut edges, invalid ring lines). Each are a geometry collection.Dangles are edges which have one or both ends which are not incident on another edge endpoint. Cut edges are connected at both ends but do not form part of polygon. Invalid ring lines form rings which are invalid (bowties, etc).

Once you have individual polygons, you can use polygon.area() to calculate area

Edit 1:

Based on your comment, I am providing a sample code with a is similar set of polylines. It performs following steps

  1. takes line1 coords and line2 coords as input
  2. Creates a combined array of coordinates from line 1 and reversed coordinates of line 2
  3. Creates a linestring out of combined coords from step 2
  4. Performs unary_union operation and then gets multilinestrings as output
  5. Passes Multilinestring obtained from step 4 to polygonize() method
  6. Polygonize returns a list of polygons which are then converted to EPSG:3857 for calculating their area in meters.
  7. Saves the output geojson to file
from shapely.ops import polygonize, unary_union
import numpy as np
from shapely.geometry import LineString
import geopandas
import pprint

line1_coords = [[72.83068656921387,18.93860101523069],[72.82768249511719,18.94176717714768],
                [72.82570838928223, 18.944770916282575],[72.82536506652832,18.94753106134636],
                [72.82759666442871,18.950778232391926],[72.83034324645996,18.953944163322795],
                [ 72.83317565917969, 18.955324165697053],[72.83532142639159,18.95824648598139],
                [72.83532142639159,18.962386352071892],[72.83317565917969,18.96490269104391],
                [72.83025741577148,18.967175480699648],[ 72.82716751098633,18.97009759329888],
                [72.82682418823242, 18.97261381587748]]

line2_coords =[[72.82716751098633,18.938195092694194],[72.83068656921387,18.945501547359086],
               [72.8324031829834,18.949641729713196],[72.83145904541016,18.95256414953349],
               [72.83068656921387,18.956298278146946],[72.8293991088867,18.96181814124214],
               [72.83008575439453,18.96652611538819],[72.83017158508301,18.968393033841075],
               [72.83060073852539,18.972694983715265],[72.83042907714844,18.97383132929355]]

#lines = [line1_coords, line2_coords]

combined_coords = line1_coords
line2_coords.reverse()
combined_coords.extend(line2_coords)
combined_coords.append(combined_coords[0])

combined_polyline = LineString(combined_coords)
mls = unary_union(combined_polyline)

polygons = list(polygonize(mls))
gdf = geopandas.GeoDataFrame(geometry=polygons, crs=4326)
gdf = gdf.to_crs(3857)
print(gdf.crs)


total_area = 0

for i, polygon in enumerate(gdf.geometry.to_list()):
    total_area = total_area + polygon.area
    print("Area of polygon# " + str(i)+ ": " + str(polygon.area) + "sq mtrs")

print("Total Area: " + str(total_area) +" sq mtrs")


gdf.to_file(open("out_geom.json", "wb"), driver='GeoJSON')

The above script produces following output. Check screenshot for polygons identified.

epsg:3857
Area of polygon# 0: 59422.22009810131sq mtrs
Area of polygon# 1: 663689.3374683487sq mtrs
Area of polygon# 2: 659717.8415444194sq mtrs
Area of polygon# 3: 202110.5252678189sq mtrs
Total Area: 1584939.9243786884 sq mtrs

enter image description here

5
  • I replaced the lines with: lines = [ ((0, 0), (1, 1),(2,0)), ((0, 0), (0, -1), (2,1)), ((0, 1), (1, -1)) ], it won't work. Not sure what is the requirements of this function.
    – Chang
    Jan 16, 2022 at 4:08
  • provided a sample working script for reference
    – Abhilshit
    Jan 16, 2022 at 23:57
  • Thanks. The solution is very close to what I needed. I realized the flowline must share a vertex. In your case, you close the opening at both ends. In my case, I want to remove the head and tail. Second, it seems does not work when there are multiple flowlines are together., so one flowline is shared by several polygons.
    – Chang
    Jan 22, 2022 at 1:27
  • I found a commercial tool that actually solves this: ian-ko.com/ET_GeoWizards/ETGW12UG_WEB/etgw/convertPl2Pg.htm
    – Chang
    Jan 22, 2022 at 17:30
  • A similar solution was provided here: rocksandwater.net/blog/2018/08/polygons-from-lines
    – Chang
    Jan 24, 2022 at 2:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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