2

I am working on a shapefile in python using geopandas and gdal. I am looking to create meshgrid (with regular 1000m interval points) inside the polygon shapefile. I have reprojected the file so that units can be meters. However, I could not find any direct way to implement this. Can any one guide in this regard?

I am sharing the code, I have tried so far:

from osgeo import gdal, ogr
import numpy as np
import matplotlib.pyplot as plt
import os
import sys
import pandas as pd
import geopandas as gpd
from shapely.geometry import Polygon
source_ds = ogr.Open(r"E:\review paper\sample tb data for recon\descend\tiffbt\alaska_bound.shp")
boundFile =gpd.read_file(r"E:\review paper\sample tb data for recon\descend\tiffbt\alaska_bound.shp")
bound_project = boundFile.to_crs({'init': 'EPSG:3572'})
print(bound_project.crs)
print(bound_project.total_bounds)

The coordinate system and bounding box coordinates are as below (output of above code):

+init=epsg:3572 +type=crs
[-2477342.73003557 -3852592.48050272  1305143.81797914 -2054961.64359753]
0

1 Answer 1

2

Find xmin, xmax, ymin, ymax. Create ranges from xmin to xmax and ymin to ymax with the interval of your spacing. Combine all possible coordinates using meshgrid:

import numpy as np
import geopandas as gpd

shapefile = r"C:\GIS\data\testdata\ak_riks.shp"
spacing = 100000

df = gpd.read_file(shapefile)
xmin, ymin, xmax, ymax = df.total_bounds

xcoords = [i for i in np.arange(xmin, xmax, spacing)]
ycoords = [i for i in np.arange(ymin, ymax, spacing)]

pointcoords = np.array(np.meshgrid(xcoords, ycoords)).T.reshape(-1, 2) #A 2D array like [[x1,y1], [x1,y2], ...
points = gpd.points_from_xy(x=pointcoords[:,0], y=pointcoords[:,1])
grid = gpd.GeoSeries(points, crs=df.crs)
grid.name = 'geometry'

#If you just want to points inside polygons:
gridinside = gpd.sjoin(gpd.GeoDataFrame(grid), df[['geometry']], how="inner")

#Plot
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(15, 15))
df.plot(ax=ax, alpha=0.7, color="pink", edgecolor='red', linewidth=1)
grid.plot(ax=ax, markersize=150, color="blue")
gridinside.plot(ax=ax, markersize=100, color="yellow")

enter image description here

2
  • At the line: grininside= gpd.sjoin(gpd.GeoDataFrame(grid), df[['geometry']], how="inner") , I am getting the following error. Spatial indexes require either rtree or pygeos. See installation instructions at geopandas.org/install.html Could you please guide me with this. I have installed both rtree and pygeos, and have geopandas as well.
    – srinivas
    Oct 14, 2022 at 10:25
  • 1
    I am able to resolve it by restarting the kernel. Thanks.
    – srinivas
    Oct 14, 2022 at 12:27

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.