Script to only keep waypoints within area?

Hello,

I need to loop through a big set of waypoints that cover a whole nation.

To speed things up, I’d like to first filter the set and only keep waypoints that lie within a region, eg. “filter.py waypoints.gpx region1.gpx → region1.gpx”.

Would someone have a script handy?

Thank you.

Here’s one way:

import json
from shapely.geometry import Point, Polygon, shape
import gpxpy
import gpxpy.gpx
import glob
import os

"""
Loop through list of locations, only keep those that lie within a region (polygon), and save as GPX file

Pre-requisite: Regions as geojson files
"""

#Read locations
locations_file = open('waypoints.gpx', mode='rt', encoding='utf-8')
locations_gpx = gpxpy.parse(stations_file)
locations_file.close()
print("Number of locations: ",len(locations_gpx.waypoints))

#Loop through regions
for file in glob.glob('*.geojson'):
	print("Handling ", file)

	#Output to GPX
	BASENAME, EXTENSION = os.path.splitext(os.path.basename(file))
	f = open(f"{BASENAME}.locations.gpx", 'w',encoding='utf-8')
	outputgpx = gpxpy.gpx.GPX()

	with open(file,encoding='utf-8') as f:
		data = json.load(f)
	#Read region
	for feature in data['features']:
		polygon = shape(feature['geometry'])
		
		#Loop through locations
		for waypoint in locations_gpx.waypoints:
			location_coords = Point(waypoint.longitude,waypoint.latitude)

			if polygon.contains(location_coords):
				print(f"Location included: {waypoint.name}")
				wpt = gpxpy.gpx.GPXWaypoint(waypoint.latitude, waypoint.longitude)
				wpt.name = waypoint.name
				outputgpx.waypoints.append(wpt)

	output = outputgpx.to_xml()
	f.write(output)
	f.close()