Create point grid inside a shapefile using python

Question

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]

Answer 1

It's not clear if you are trying to create a grid of boxes or a grid of points. To change to points use:

# create a grid for geometry
gdf_grid = gpd.GeoDataFrame(
    geometry=[
        shapely.geometry.Point(x, y)
        for x in np.arange(a, c, STEP)
        for y in np.arange(b, d, STEP)
    ],
    crs=crs,
).to_crs(gdf.crs)

• have used 50km instead of 1000m for demonstration purposes
• with Alaska it for polygons it is necessary to take into account the antimeridian. Without this you will have polygons than span in excess of 350 degrees when re-projected to EPSG:4326
• approach is simple
  1. obtain Alaska geometry shape file
  2. project to a CRS in meters. Have used UTM
  3. get total_bounds
  4. construct grid of geometry objects using 3
  5. restrict grid of geometry to ones that intersect with geometry
• you will observe at such latitudes there will be distortion between UTM and EPSG:4326 as expected (the nature of projections)

full code

import geopandas as gpd
import numpy as np
import shapely.geometry

gdf = gpd.read_file("https://www2.census.gov/geo/tiger/TIGER2018/ANRC/tl_2018_02_anrc.zip")

STEP = 50000
crs = gdf.estimate_utm_crs()
# crs = "EPSG:3338"
a, b, c, d = gdf.to_crs(crs).total_bounds

# create a grid for geometry
gdf_grid = gpd.GeoDataFrame(
    geometry=[
        shapely.geometry.box(minx, miny, maxx, maxy)
        for minx, maxx in zip(np.arange(a, c, STEP), np.arange(a, c, STEP)[1:])
        for miny, maxy in zip(np.arange(b, d, STEP), np.arange(b, d, STEP)[1:])
    ],
    crs=crs,
).to_crs(gdf.crs)

# exclude geometries that cross antimeridian 
gdf_grid = gdf_grid.loc[~gdf_grid["geometry"].bounds.pipe(lambda d: d["maxx"] - d["minx"]).ge(350)]

# restrict grid to only squares that intersect with geometry
gdf_grid = (
    gdf_grid.sjoin(gdf.dissolve().loc[:,["geometry"]])
    .pipe(lambda d: d.groupby(d.index).first())
    .set_crs(gdf.crs)
    .drop(columns=["index_right"])
)
m = gdf.explore(color="red", style_kwds={"fillOpacity":0})
gdf_grid.explore(m=m)

output