In [42]:
import mmi
import numpy as np
import shapely.geometry
import requests
import zmq
import zmq.eventloop
import osgeo.osr
import shapely.geometry
In [43]:
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)
In [44]:
names = {"xk", "yk", "flowelemnode"}
for var in names:
mmi.send_array(pushstream, metadata={"get_var": var})
In [45]:
arrs = {}
for name in names:
A, metadata = mmi.recv_array(sub)
arrs[metadata["name"]] = A
In [46]:
arrs
Out[46]:
In [47]:
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
Out[47]:
In [69]:
# 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)]
In [70]:
polys = [shapely.geometry.mapping(shapely.geometry.Polygon(xy)) for xy in polycoords]
In [71]:
features = [dict(type="Feature", geometry=poly, properties=dict(index=i)) for i, poly in enumerate(polys)]
collection = dict(type="FeatureCollection", features=features)
In [72]:
import json
json.dump(collection, open("grid.json","w"))
In [68]:
Out[68]:
In [ ]: