mask_geojson.py

import json
import numpy as np
import matplotlib.pyplot as plt
from pyproj import Proj, transform

wgs84_proj = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "
au_albers = "+proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs "
# Canberra (Australia)
lat = -35.28
lon = 149.13
inProj = Proj(wgs84_proj)
outProj = Proj(au_albers)
x, y = transform(inProj,outProj,lon,lat)
print(x, y)
exit()

def is_point_in_path(x, y, poly):
    
    num = len(poly)
    i = 0
    j = num - 1
    c = False
    for i in range(num):
        if ((poly[i][1] > y) != (poly[j][1] > y)) and \
                (x < poly[i][0] + (poly[j][0] - poly[i][0]) * (y - poly[i][1]) /
                                  (poly[j][1] - poly[i][1])):
            c = not c
        j = i
    return c


poly = """{
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [100.0, 10.0], [101.0, 10.0], [101.5, 11.0],
            [100.0, 11.0], [100.0, 10.0]
          ]
        ]
      },
      "properties": {
        "prop0": "value0",
        "prop1": { "this": "that" }
      }
}"""

print(poly)
p = json.loads(poly)
coords = p['geometry']['coordinates'][0]

print(is_point_in_path(100.5, 10.5, coords))
print(is_point_in_path(101.5, 9.5, coords))

canvas = np.zeros((300, 300))

for j in range(300):
    b = 9. + (j/100.)
    for i in range(300):
        a = 99. + (i/100.)
        if is_point_in_path(a, b, coords):
            print("write", j,i)
            canvas[j,i] = 1

plt.imsave("out.png", canvas)