codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
#,!/usr/bin/env python import math from functools import partial from operator import itemgetter import fiona import pyproj from matplotlib import pyplot from ujson import dumps # noqa from shapely.ops import transform from descartes import PolygonPatch from shapely.geometry import mapping, shape # noqa import shapely.speedups shapely.speedups.available and shapely.speedups.enable() def meters_per_pixel(zoom, lat): """ ground resolution = cos(latitude * pi/180) * earth circumference / map width """ return (math.cos(lat * math.pi/180.0) * 2 * math.pi * 6378137) / (256 * 2**zoom) def map_scale(zoom, lat, dpi=96.0): """ map scale = 1 : ground resolution * screen dpi / 0.0254 meters/inch """ res = meters_per_pixel(zoom, lat) return (res * dpi) / 0.0254 def zoom_to_scaledenom(z1, z2=False): """ Converts zoom level to mapnik's scaledenominator pair for EPSG:3857 """ if not z2: z2 = z1 s = 279541132.014 z1 = (s / (2 ** (z1 - 1)) + s / (2 ** (z1 - 2))) / 2 z2 = (s / (2 ** (z2 - 1)) + s / (2 ** z2)) / 2 # return 100000000000000, 1 return z1, z2 def pixel_size_at_zoom(z, l=1): """ Converts l pixels on tiles into length on zoom z """ return int(math.ceil(l * 20037508.342789244 / 256 * 2 / (2 ** z))) zoom = 3 lat = -28.226970038918328 mpp = meters_per_pixel(zoom, lat) scale = map_scale(zoom, lat) ds = fiona.open("data/regions/LGA", "r") fs = list(ds) gs = map(shape, filter(None, map(itemgetter("geometry"), fs))) epsg_4283_to_epsg_3857 = partial( pyproj.transform, pyproj.Proj(init="epsg:4283"), pyproj.Proj(init="epsg:3857") ) epsg_3857_to_epsg_4326 = partial( pyproj.transform, pyproj.Proj(init="epsg:3857"), pyproj.Proj(init="epsg:4326") ) g1 = gs[510] # largest geometry in LGA regions ~165KB GeoJSON g2 = transform(epsg_4283_to_epsg_3857, g1) #g3 = g2.simplify(pixel_size_at_zoom(zoom, 0.5)) g3 = g2.buffer(pixel_size_at_zoom(zoom, 10)).simplify(pixel_size_at_zoom(zoom, 3)).buffer(pixel_size_at_zoom(zoom, -10)) g4 = transform(epsg_3857_to_epsg_4326, g3) print "size(g1):", len(dumps(mapping(g1))) print "size(g4):", len(dumps(mapping(g4))) #Map resolution = 156543.04 meters/pixel * cos(latitude) / (2 ^ zoomlevel) #Komzzzpa> well, algo is: # 1) convert to EPSG:3857 # 2) make simplification for 0.5px with .simplify(pixel_size_at_zoom(0.5)) # 3) convert to EPSG:4326 # 4) serialize/show #<Komzzzpa> for polygons, buffer(10px).simplify(3px).buffer(-10px) might be a lot better #<Komzzzpa> for adjacent polygons, whole other ideas are needed :3 # Create a matplotlib figure fig = pyplot.figure(num=1, figsize=(10, 4), dpi=180) # Create a subplot ax = fig.add_subplot(121) polys = g4.geoms if g4.type == "MultiPolygon" else [g4] for polygon in polys: # Make the polygon into a patch and add it to the subplot patch = PolygonPatch(polygon, facecolor='#cccccc', edgecolor='#999999') ax.add_patch(patch) # Fit the figure around the polygon's bounds, render, and save minx, miny, maxx, maxy = polygon.bounds w, h = maxx - minx, maxy - miny ax.set_xlim(minx - 0.2*w, maxx + 0.2*w) ax.set_ylim(miny - 0.2*h, maxy + 0.2*h) ax.set_aspect(1) pyplot.show() #fig.savefig('patches.png')
Private
[
?
]
Run code
Submit