import mmi
import numpy as np
import shapely.geometry
import requests
import zmq
import zmq.eventloop
import osgeo.osr
import shapely.geometry
ctx = zmq.Context()
models = requests.get('http://localhost:22222/models').json()
model = models[0]
node = "localhost"
pubport = model["ports"]["PUB"]
sub = ctx.socket(zmq.SUB)
addr = "tcp://{node:s}:{port:d}".format(node=node, port=pubport)
sub.connect(addr)
sub.setsockopt_string(zmq.SUBSCRIBE, u'')
pullport = model["ports"]["PULL"]
push = ctx.socket(zmq.PUSH)
addr = "tcp://{node:s}:{port:d}".format(node=node, port=pullport)
push.connect(addr)
pushstream = zmq.eventloop.zmqstream.ZMQStream(push)
names = {"xk", "yk", "flowelemnode"}
for var in names:
mmi.send_array(pushstream, metadata={"get_var": var})
arrs = {}
for name in names:
A, metadata = mmi.recv_array(sub)
arrs[metadata["name"]] = A
arrs
{'flowelemnode': array([[ 780, 853, 781, -1, -1, -1], [ 781, 853, 13112, -1, -1, -1], [ 781, 13112, 13110, -1, -1, -1], ..., [28792, 28797, 28796, 28798, -1, -1], [28793, 28798, 28796, 28794, -1, -1], [28780, 28783, 28784, 28786, 28782, 28781]], dtype=int32), 'xk': array([ 547994.85513183, 547854.78975898, 547712.6433233 , ..., 595202.07393365, 595315.15648816, 595225.4700929 ]), 'yk': array([ 4204745.95785394, 4204848.64553312, 4204954.42679592, ..., 4214406.36377765, 4214530.90081606, 4214617.57328289])}
flowelemnode = np.ma.masked_equal(arrs["flowelemnode"], -1)
xk = arrs["xk"]
yk = arrs["yk"]
points = np.c_[xk, yk]
src_srs = osgeo.osr.SpatialReference()
src_srs.ImportFromEPSG(model["epsg"])
dst_srs = osgeo.osr.SpatialReference()
dst_srs.ImportFromEPSG(4326)
transform = osgeo.osr.CoordinateTransformation(src_srs, dst_srs)
wkt_points = np.array(transform.TransformPoints(points))
wkt_points
array([[-122.45342788, 37.98909466, 0. ], [-122.45501605, 37.99002753, 0. ], [-122.45662776, 37.99098838, 0. ], ..., [-121.91460836, 38.0724204 , 0. ], [-121.91330272, 38.07353071, 0. ], [-121.91431353, 38.07432118, 0. ]])
# compute number of nodes per element
elemtype = (~flowelemnode.mask).sum(1)
# compute the places where we can split the 1d array
splitidx = np.cumsum(np.r_[elemtype][:-1])
# flatten the node numbers, remove the missings
idx = flowelemnode[(~flowelemnode.mask)]-1
# use the indices to lookup lon, lat
# only keep x and y, transform to list
polycoords = [x[:,:2] for x in np.split(wkt_points[idx],splitidx)]
polys = [shapely.geometry.mapping(shapely.geometry.Polygon(xy)) for xy in polycoords]
features = [dict(type="Feature", geometry=poly, properties=dict(index=i)) for i, poly in enumerate(polys)]
collection = dict(type="FeatureCollection", features=features)
import json
json.dump(collection, open("grid.json","w"))
[array([[-122.4806033 , 38.09635956], [-122.47300867, 38.10211462], [-122.48402421, 38.09970107]]), array([[-122.48402421, 38.09970107], [-122.47300867, 38.10211462], [-122.48052373, 38.10561733]]), array([[-122.48402421, 38.09970107], [-122.48052373, 38.10561733], [-122.48459261, 38.10532472]]), array([[-122.48402421, 38.09970107], [-122.48459261, 38.10532472], [-122.48897247, 38.1033987 ]]), array([[-122.48897247, 38.1033987 ], [-122.48459261, 38.10532472], [-122.48771973, 38.10716098]]), array([[-122.47300867, 38.10211462], [-122.46462764, 38.1077655 ], [-122.4757648 , 38.10873706]]), array([[-122.47300867, 38.10211462], [-122.4757648 , 38.10873706], [-122.48052373, 38.10561733]]), array([[-122.46462764, 38.1077655 ], [-122.47193063, 38.11491708], [-122.4757648 , 38.10873706]]), array([[-122.46462764, 38.1077655 ], [-122.45617388, 38.11251683], [-122.45986205, 38.11667626]]), array([[-122.46462764, 38.1077655 ], [-122.45986205, 38.11667626], [-122.47193063, 38.11491708]])]