[MapProxy] Re: Using hex tile indices (ArcGIS tile cache)

[Josh Doe]

I've made a little Python script to assist in generating a MapProxy
grid from the JSON response of an ArcGIS server.

Since the origin of MapProxy grids are lower-left (sw), and ArcGIS
uses upper-left (nw) as its origin, you have to calculate the other
corner (well, at least the bottom). Also, since some resolution sets
aren't neat multiples of 2, you have to find the least common multiple
of the resolutions, and then find the multiple of them which at least
includes the extent. It's not hard, but it took me a little while to
get it right, so hopefully others can benefit from this.

I've attached the script, and you can run it like:

wget -O mapserver_json
http://gismaps.virginia.gov/arcgis2/rest/services/VBMP2006_2007/MapServer?f=json
python arcgisjson2mapproxygrid.py mapserver_json

Use the generated grid in your mapproxy.yaml, and be sure your source
uses the grid and specifies 'origin: nw'.

Regards,
-Josh

On Mon, Mar 21, 2011 at 3:52 AM, Oliver Tonnhofer <olt at omniscale.de> wrote:
>
> On 19.03.2011, at 06:03, Josh Doe wrote:
>
>> I've added this to my fork [1].
>
> Great. Thanks for the small but complete patch. I pulled your changes into the default branch.
>
>
> Regards,
> Oliver
>
> --
> Oliver Tonnhofer    | Omniscale GmbH & Co KG    | http://omniscale.de
> http://mapproxy.org | https://bitbucket.org/olt
>
>
>
>
-------------- next part --------------
import json
import math
import string
import sys

# from J.F. Sebastian: http://stackoverflow.com/users/4279/j-f-sebastian
# from this question: http://stackoverflow.com/questions/147515
def gcd(a, b):
    """Return greatest common divisor using Euclid's Algorithm."""
    while b: a, b = b, a % b
    return a
def lcm(a, b):
    """Return lowest common multiple."""
    return a * b // gcd(a, b)
def lcmm(args):
    """Return lcm of args."""
    return reduce(lcm, args)


if len(sys.argv) < 2:
    print "Usage: %s MAPSERVER_JSON" % (sys.argv[0])
    sys.exit(0)

# Read MapServer service description via json
f = open(sys.argv[1], 'r')
d = json.load(f)
f.close()

# get extent
fe = d['fullExtent']
xmin = fe['xmin']
ymin = fe['ymin']
xmax = fe['xmax']
ymax = fe['ymax']
fe_srs = fe['spatialReference']

# get basic tile info
ti = d['tileInfo']
tw = ti['rows']
th = ti['cols']
fmt = ti['format']
left = ox = ti['origin']['x']
top = oy = ti['origin']['y']
ti_srs = ti['spatialReference']

if ti_srs != fe_srs:
  print 'Error: SRS of fullExtent does not match SRS of tileInfo'
  sys.exit(-1)

# get resolutions and scales
a=[]
scale=[]
for l in d['tileInfo']['lods']:
  a.append(str(l['resolution']))
  scale.append(l['scale'])
res = string.join(a,',')

# find least common multiple of all scales and convert to resolution
inc = lcmm(scale)/scale[0]*d['tileInfo']['lods'][0]['resolution']

# find multiple of inc which encompasses extent
mult = math.ceil((oy-ymin)/(inc*th))
bottom = oy-mult*inc*th
mult = math.ceil((xmax-ox)/(inc*tw))
right = ox+mult*inc*tw

print """
#For the source use one of the following URLs, changing basepath and format:
#  http://arcgisserver/MapServer/tile/%(z)s/%(y)s/%(x)s
#  http://arcgiscache/%(arcgiscache_path)s.jpg
#
#Also be sure to use 'origin: nw'"""

print """
grids:
  arcgiscache_grid:
    srs: #FIXME: find SRS (WKID) that corresponds to WKT
    tile_size: [%d, %d]
    res: [%s]
    bbox: [%f, %f, %f, %f]""" % (tw, th, res, left, bottom, right, top)