1

I have two geospatial shapefiles.

1. The first one is a contour shapefile, generated via following code

rasterDs = gdal.Open('raster.tif')
rasterBand = rasterDs.GetRasterBand(1)
contourPath="contour.shp"
contourDs = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(contourPath)

#define layer name and spatial 
contourShp = contourDs.CreateLayer('contour', proj)

#define fields of id and elev
fieldDef = ogr.FieldDefn("ID", ogr.OFTInteger)
contourShp.CreateField(fieldDef)
fieldDef = ogr.FieldDefn("elev", ogr.OFTReal)
contourShp.CreateField(fieldDef)

#The raster band is taken from the first cell
gdal.ContourGenerate(rasterBand, 12.5, 240, [], 1, 0., 
                     contourShp, 0, 1)

contourDs.Destroy()

2. The second one is a path file, containing LineString and have been generated by following code

import fiona
import pandas as pd

def generate_shp_file(line_csv, des_output):
    """
        Converting the final coordinates stored as csv to shapefile.
    """
    #import dataframe from csv file
    coordinate_df = pd.read_csv(line_csv, header = 0)

    # define schema
    schema = {
                'geometry':'LineString',
                'properties':[('Name', 'str')]
            }
    
    lineShp = fiona.open(des_output, mode = 'w', driver = 'ESRI Shapefile',
                         schema = schema, crs = "EPSG:4326")

    # Get the list of points
    coordinates = []
    rowName = ''
    for index, row in coordinate_df.iterrows():
        coordinates.append((row.Y,row.X))
        rowName = row.Name

    # Save record according to the subject of shapefile.
    rowDict = {
                    'geometry': {
                                    'type':'LineString',
                                    'coordinates': coordinates
                                },
                    'properties': {
                                    'Name' : rowName
                                },
              }
    lineShp.write(rowDict)
    lineShp.close()

line_csv file is a csv file which contains latitude and longitude under header X and Y respectively.

The task is to find the intersection points between the contour shapefile and path shapefile.

I am using following code to do so

import geopandas as gpd
import warnings
# Load contour and linestring shapefiles
contour = gpd.read_file(contourPath)
linestring = gpd.read_file(shapefile_path)

print(contour.crs)
# # Reproject both geometries to EPSG:4326
contour = contour.to_crs(epsg=32645)
linestring = linestring.to_crs(epsg=32645)

intersections = linestring.intersection(contour)
# Check if any intersections were found
if not intersections.empty:
    print("Intersections found between linestring and contour.")
    intersections.to_file('intersections.shp', crs=contour.crs)
else:
    print("No intersections found. There might be gaps or misalignments.")

The code is not able to find intersections and returns empty LINESTRING. However, it is quite evident in the below image that both these vectors are intersecting at least at 3 points. Also, I am getting a warning geopandas\base.py:31: UserWarning: The indices of the two GeoSeries are different.

Why is the intersection list is empty?

Extract one point from the LineString geometry`

Picture showing intersection of contour and line at 3 points

2
  • If you are sure that the code is correct and the result us null can you make sure that if the CRS is the same?
    – GforGIS
    Sep 12, 2023 at 21:22
  • In this line I have made the CRS same contour = contour.to_crs(epsg=32645) linestring = linestring.to_crs(epsg=32645) Sep 13, 2023 at 3:59

1 Answer 1

2

Use gpd.overlay with keep_geom_type=False:

import geopandas as gpd

contour = gpd.read_file(r"C:\GIS\data\testdata\contour.shp")
river = gpd.read_file(r"C:\GIS\data\testdata\rivers.shp")
river = river.to_crs(contour.crs)

#i = contour.intersection(river)
# i.unique()
# <GeometryArray>
# [<LINESTRING Z EMPTY>]

i2 = gpd.overlay(df1=contour, df2=river, how="intersection", keep_geom_type=False)
# i2.geometry.head(2)
# 0    POINT (431730.187 6370770.016)
# 1    POINT (431227.744 6370352.271)

ax = contour.plot(linewidth=1, edgecolor="gray", figsize=(20,20))
river.plot(ax=ax, linewidth=1.5, edgecolor="blue")
i2.plot(ax=ax, markersize=100, marker="o", color="red")

enter image description here

2
  • 1
    Thanks it worked, can you please explain or guide me to some tutorial about it. Thanks Sep 13, 2023 at 9:14
  • 1
    Nice! See: Set operations with overlay. Intersecting two line layers will result in points if the lines doesnt align exactly.
    – BERA
    Sep 13, 2023 at 11:00

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.