Last active
January 5, 2023 18:09
-
-
Save dhardy92/28bba665235b4d7e06e0445b22778f64 to your computer and use it in GitHub Desktop.
Generate a Squadrat data GeoJSON format file
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 | |
}) | |
}); |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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)) |
New version add ubersquadrat and ubersquadratinhos based on https://www.baeldung.com/cs/max-size-square-matrix#two-pointers-approach
modified main.js file added
this file is created in directory òpenlayerin the Openlayer projet with command
npm create ol-app openlayer`. simply replace it.
add clusters and yard* display approching squadrats style inspired by https://www.geeksforgeeks.org/find-the-number-of-islands-using-dfs/
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
The file is a GeoJSON format that embede 2 MultiPolygon features :
This file can now be used as source in inrerface build with OpenLayers for example