Skip to content

Instantly share code, notes, and snippets.

@dhardy92

dhardy92/main.js

Last active January 5, 2023 18:09
Show Gist options
  • Save dhardy92/28bba665235b4d7e06e0445b22778f64 to your computer and use it in GitHub Desktop.
Save dhardy92/28bba665235b4d7e06e0445b22778f64 to your computer and use it in GitHub Desktop.
Download ZIP
Generate a Squadrat data GeoJSON format file
Raw
main.js
import './style.css';
import {Map, View} from 'ol';
import TileLayer from 'ol/layer/WebGLTile';
import OSM from 'ol/source/OSM';
import {useGeographic} from 'ol/proj.js';
import {Vector as VectorLayer} from 'ol/layer.js';
import VectorSource from 'ol/source/Vector.js';
import {Fill, Stroke, Style} from 'ol/style.js';
import GeoJSON from 'ol/format/GeoJSON.js';
const styles = {
'ubersquadrat': new Style({
stroke: new Stroke({
color: 'red',
width: 2,
}),
}),
'ubersquadratinho': new Style({
stroke: new Stroke({
color: 'brown',
width: 2,
}),
}),
'squadratinhos': new Style({
fill: new Fill({
color: 'rgba(255, 200, 180, 0.5)',
}),
}),
'squadrats': new Style({
fill: new Fill({
color: 'rgba(255, 200, 180, 0.5)',
}),
}),
'yardinho': new Style({
fill: new Fill({
color: 'rgba(255, 135, 55, 0.5)',
}),
}),
'yard': new Style({
fill: new Fill({
color: 'rgba(255, 255, 255, 0.5)',
}),
}),
};
useGeographic();
const size = 256;
const canvas = document.createElement('canvas');
canvas.width = size;
canvas.height = size;
const context = canvas.getContext('2d');
context.strokeStyle = 'white';
context.textAlign = 'center';
context.font = '24px sans-serif';
const lineHeight = 30;
const vector = new VectorLayer({
source: new VectorSource({
url: 'geojson/squadrats.json',
format: new GeoJSON(),
}),
style: function (feature) {
return styles[feature.getId()];
},
});
const map = new Map({
target: 'map',
layers: [
new TileLayer({
source: new OSM()
}),
vector,
],
view: new View({ //paris
center: [2.3488, 48.85341],
zoom: 14
})
});
Raw
squadrats.py
#usage: squadrats.py [-h] [-f FILE]
#
#options:
# -h, --help show this help message and exit
# -f FILE, --file FILE GPX input file or directory containing GPX files
import gpxpy
import gpxpy.gpx
import math
import sys
import argparse
import json
import os
#augment global variable set tiles{} with new tiles for the current GPX file
def dofile(file):
gpx_file = open(file, 'r')
gpx = gpxpy.parse(gpx_file)
global bounds
b = gpx.get_bounds()
if bounds != []:
bounds[0] = max (bounds[0],b.max_latitude) #north
bounds[1] = min (bounds[1],b.min_latitude) #south
bounds[2] = max (bounds[2],b.max_longitude) #east
bounds[3] = min (bounds[3],b.min_longitude) #west
else:
bounds = [b.max_latitude,b.min_latitude,b.max_longitude,b.min_longitude]
for track in gpx.tracks:
for segment in track.segments:
for point in segment.points:
for zoom in (14,17):
tile = deg2num(point.latitude, point.longitude, zoom)
tiles.add('{0}/{1}/{2}'.format(zoom, tile[0], tile[1]))
#calculate coordinate of the tile for the current point of the current GPX file
def deg2num(lat_deg, lon_deg, zoom):
lat_rad = math.radians(lat_deg)
n = 2.0 ** zoom
xtile = int((lon_deg + 180.0) / 360.0 * n)
ytile = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
return (xtile, ytile)
#get coord from tile number (shape)
def num2deg(xtile, ytile, zoom):
n = 2.0 ** zoom
lon_deg = xtile / n * 360.0 - 180.0
lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
lat_deg = math.degrees(lat_rad)
return (lon_deg, lat_deg)
def getsquare(xtile, ytile, zoom, size=1):
return [ num2deg(xtile, ytile, zoom),
num2deg(xtile, ytile+size, zoom),
num2deg(xtile+size, ytile+size, zoom),
num2deg(xtile+size, ytile, zoom),
num2deg(xtile, ytile, zoom) ]
#function determine if a submatrix (rectangle) is full.
# based on https://www.baeldung.com/cs/max-size-square-matrix
def isFull(tl,br,z):
global prefix
area = ( br['x'] - tl['x'] + 1) * (br['y'] - tl['y'] + 1)
s = prefix[z][br['y']][br['x']]
s = s - prefix[z][tl['y']-1][br['x']]
s = s - prefix[z][br['y']][tl['x']-1]
s = s + prefix[z][tl['y']-1][tl['x']-1]
return s == area
# Program to count islands in boolean 2D matrix
class Graph:
def __init__(self, g):
self.ROW = len(g)
self.COL = len(g[0])
self.graph = g
# A function to check if a given cell
# (row, col) can be included in DFS
def isSafe(self, i, j, visited):
# row number is in range, column number
# is in range and value is 1
# and not yet visited
return (i >= 0 and i < self.ROW and
j >= 0 and j < self.COL and
not visited[i][j] and self.graph[i][j])
# A utility function to do DFS for a 2D
# boolean matrix. It only considers
# the 8 neighbours as adjacent vertices
def DFS(self, i, j, visited):
# These arrays are used to get row and
# column numbers of 4 cardinal neighbours
# of a given cell
rowNbr = [-1, 0, 0, 1]
colNbr = [ 0, -1, 1, 0]
# Mark this cell as visited
visited[i][j] = True
c=[(j,i)]
# Recur for all connected neighbours
for k in range(4):
if self.isSafe(i + rowNbr[k], j + colNbr[k], visited):
c = c + self.DFS(i + rowNbr[k], j + colNbr[k], visited)
return c
# The main function that returns
# count of islands in a given boolean
# 2D matrix
def getIslands(self):
# Make a bool array to mark visited cells.
# Initially all cells are unvisited
visited = [[False for j in range(self.COL)]for i in range(self.ROW)]
# Initialize count as 0 and traverse
# through the all cells of
# given matrix
count = 0
clusters = []
for i in range(self.ROW):
for j in range(self.COL):
# If a cell with value 1 is not visited yet,
# then new island found
if visited[i][j] == False and self.graph[i][j] == 1:
# Visit all cells in this island
# and increment island count
clusters.append(self.DFS(i, j, visited))
count += 1
return clusters
#main
parser = argparse.ArgumentParser()
parser.add_argument("-f", "--file", help="GPX input file or directory containing GPX files")
args = parser.parse_args()
tiles = set()
bounds = [] #[north,south,east,west]
if os.path.isdir(args.file):
for file in os.scandir(args.file):
if file.path.endswith('gpx'):
dofile(file)
else:
dofile(args.file)
geojson = {
"type": "FeatureCollection",
"features":[
{
"type": "Feature",
"id": "squadratinhos",
"title": "Squadratinhos",
"geometry": {
"type": "MultiPolygon",
"coordinates":[],
},
"properties": {}
},
{
"type": "Feature",
"id": "squadrats",
"title": "Squadrats",
"geometry": {
"type": "MultiPolygon",
"coordinates":[],
},
"properties": {}
},
{
"type": "Feature",
"id": "ubersquadratinho",
"title": "Übersquadratinho",
"geometry": {
"type": "Polygon",
"coordinates":[],
},
"properties": {}
},
{
"type": "Feature",
"id": "ubersquadrat",
"title": "Übersquadrat",
"geometry": {
"type": "Polygon",
"coordinates":[],
},
"properties": {}
},
{
"type": "Feature",
"id": "yardinho",
"title": "Yardinho",
"geometry": {
"type": "MultiPolygon",
"coordinates":[],
},
"properties": {}
},
{
"type": "Feature",
"id": "yard",
"title": "Yard",
"geometry": {
"type": "MultiPolygon",
"coordinates":[],
},
"properties": {}
},
]
}
# init matrixes and variables for max square search
bounds_tiles = {}
matrix = {}
prefix = {}
cluster= {}
square = {}
yards = {}
for z in (14,17):
bounds_tiles[z] = [ deg2num(bounds[0],bounds[3] , z), deg2num(bounds[1], bounds[2], z) ]
prefix[z] = [[0 for column in range(bounds_tiles[z][1][0] - bounds_tiles[z][0][0] + 1)] for row in range(bounds_tiles[z][1][1] - bounds_tiles[z][0][1] + 1)]
matrix[z] = [[0 for column in range(bounds_tiles[z][1][0] - bounds_tiles[z][0][0] + 1)] for row in range(bounds_tiles[z][1][1] - bounds_tiles[z][0][1] + 1)]
cluster[z] = [[0 for column in range(bounds_tiles[z][1][0] - bounds_tiles[z][0][0] + 1)] for row in range(bounds_tiles[z][1][1] - bounds_tiles[z][0][1] + 1)]
# loop on tiles to creates matrixes and geojson documents
for i in tiles:
(z,x,y) = i.split('/')
mx=int(x) - bounds_tiles[int(z)][0][0]
my=int(y) - bounds_tiles[int(z)][0][1]
matrix[int(z)][my][mx] = 1
if (z=='17'): #squadratinhos
geojson["features"][0]["geometry"]["coordinates"].append([getsquare(int(x),int(y),int(z))])
elif (z=='14'): # squadrats
geojson["features"][1]["geometry"]["coordinates"].append([getsquare(int(x),int(y),int(z))])
# search for max square in both zoom
for z in (14,17):
square[z]={"l":0,"x":0,"y":0}
for i in range(len(matrix[z])):
for j in range(len(matrix[z][0])):
#clusters[i][j] 1 if cardinal points (NSEW) are 1
if (i > 0 and i < len(matrix[z])-1 and
j > 0 and j < len(matrix[z][0])-1 and
matrix[z][i+1][j] +
matrix[z][i][j+1] +
matrix[z][i-1][j] +
matrix[z][i][j-1] == 4):
cluster[z][i][j] = 1
#prefix[i][j] is the sum of values on the square (matrix[0][0],matrix[i],[j])
if (i == j == 0) :
prefix[z][i][j] = matrix[z][i][j]
elif (j==0) :
prefix[z][i][j] = matrix[z][i][j] + prefix[z][i-1][j]
else :
prefix[z][i][j] = matrix[z][i][j] + prefix[z][i-1][j] + prefix[z][i][j-1] - prefix[z][i-1][j-1]
# based on https://www.baeldung.com/cs/max-size-square-matrix#two-pointers-approach
for i in range(len(matrix[z])):
for j in range(len(matrix[z][i])):
while (i + square[z]['l'] < len(matrix[z]) and j + square[z]['l'] < len(matrix[z][i]) and isFull( {"x":j,"y":i}, {"x":j+square[z]['l'],"y":i+square[z]['l']}, z)):
square[z] = {"l": square[z]['l']+1, "x":j, "y":i}
# reset origin from bounding
square[z]['x']=square[z]['x'] + bounds_tiles[z][0][0]
square[z]['y']=square[z]['y'] + bounds_tiles[z][0][1]
g = Graph(cluster[z])
yards[z]= g.getIslands()
yards[z].sort(key=len,reverse=True)
if (yards[z] != None and len(yards[z]) >0):
for yards_tiles in yards[z][0]:
yard = 4 if (z == 17) else 5
t = getsquare(yards_tiles[0] + bounds_tiles[z][0][0], yards_tiles[1] + bounds_tiles[z][0][1],z)
geojson["features"][yard]["geometry"]["coordinates"].append([t])
geojson["features"][yard-4]["geometry"]["coordinates"].remove([t])
# add to geojson document
geojson["features"][2]["geometry"]["coordinates"].append(getsquare(square[17]['x'],square[17]['y'],17,square[17]['l']))
geojson["features"][3]["geometry"]["coordinates"].append(getsquare(square[14]['x'],square[14]['y'],14,square[14]['l']))
print(json.dumps(geojson))
@dhardy92
Copy link
Author

The file is a GeoJSON format that embede 2 MultiPolygon features :

  • squadratinhos (all zoom 17 visited tiles)
  • squdrats (all zoom 14 visited tiles)
    This file can now be used as source in inrerface build with OpenLayers for example

@dhardy92
Copy link
Author

New version add ubersquadrat and ubersquadratinhos based on https://www.baeldung.com/cs/max-size-square-matrix#two-pointers-approach
image

@dhardy92
Copy link
Author

modified main.js file added
this file is created in directory òpenlayerin the Openlayer projet with commandnpm create ol-app openlayer`. simply replace it.

@dhardy92
Copy link
Author

dhardy92 commented Jan 5, 2023 •
edited
Loading

add clusters and yard* display approching squadrats style inspired by https://www.geeksforgeeks.org/find-the-number-of-islands-using-dfs/
image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment