-
Notifications
You must be signed in to change notification settings - Fork 0
/
bbox_text.py
97 lines (86 loc) · 3.6 KB
/
bbox_text.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
"""
This module provides functions for geographic and projected coordinate conversion,
and for checking if a point is within a shapefile's geometry.
"""
import shapefile
from shapely.geometry import shape, Point
from pyproj import Transformer
def convert_to_geographic(proj_coords, from_epsg, to_epsg=4326):
"""
Convert projected coordinates to geographic coordinates (longitude, latitude).
Parameters:
proj_coords : tuple
The projected coordinates (x, y).
from_epsg : int
The EPSG code of the projected coordinate system.
to_epsg : int, optional
The EPSG code of the geographic coordinate system (default is 4326, WGS84).
Returns:
tuple : The geographic coordinates (longitude, latitude).
"""
# Use Transformer to convert coordinates
transformer = Transformer.from_crs(f"epsg:{from_epsg}", f"epsg:{to_epsg}", always_xy=True)
lon, lat = transformer.transform(proj_coords[0], proj_coords[1])
return lon, lat
def convert_to_projected(geo_coords, from_epsg=4326, to_epsg=3857):
"""
Convert geographic coordinates (longitude, latitude) to projected coordinates.
Parameters:
geo_coords : tuple
The geographic coordinates (longitude, latitude).
from_epsg : int
The EPSG code of the geographic coordinate system (default is 4326, WGS84).
to_epsg : int, optional
The EPSG code of the projected coordinate system (default is 3857, Web Mercator).
Returns:
tuple : The projected coordinates (x, y).
"""
# Use Transformer to convert coordinates
transformer = Transformer.from_crs(f"epsg:{from_epsg}", f"epsg:{to_epsg}", always_xy=True)
x, y = transformer.transform(geo_coords[0], geo_coords[1])
return x, y
def is_point_in_shapefile(shapefile_path, point, state):
"""
Check if a point is within any of the shapes in the shapefile using spatial filtering,
and filter by state.
Parameters:
shapefile_path : str
Path to the shapefile (.shp).
point : tuple of (longitude, latitude)
The geographic point to check.
state : str
The state to filter shapes by (e.g., "GA").
Returns:
bool : True if the point is within any shape's bounding box in the specified state,
False otherwise
"""
# Read the shapefile
sf = shapefile.Reader(shapefile_path)
try:
# Get the point's longitude and latitude
x, y = convert_to_projected(point)
point_proj = Point(x, y)
# Get the field names and find the index of the state field
fields = sf.fields[1:] # Skip the first deletion flag field
field_names = [field[0] for field in fields]
state_index = field_names.index("State")
# Iterate through each shape record
for shape_record in sf.iterShapeRecords():
# Check if the shape's state matches the specified state
if shape_record.record[state_index] == state:
# Perform a spatial filter on the bounding box of the shape
# bbox = shape_record.shape.bbox
polygon = shape(shape_record.shape)
if polygon.contains(point_proj):
record_data = {
field_name: str(shape_record.record[i])
.encode('utf-8', 'ignore')
.decode('utf-8')
for i, field_name in enumerate(field_names)
}
# Return the record as a JSON object
return record_data
except KeyError:
return False
# If the point is not within any shapes in the specified state
return False