#!/usr/bin/python

# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included
# in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
# DEALINGS IN THE SOFTWARE.

#
# Generate Soviet Millitary Grid
#
# Author: Cristian Balint <cristian.balint@gmail.com>
#
#  - This script requires gdal-python module to be installed
# to run it on Fedora for example type 'yum install gdal-python'
# or use latest FWtools aviable for win32 enviroment.
#


try:
    from osgeo import ogr
except ImportError:
    import ogr


alphabet = ['A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T','U','V','Z','X']

# output files
file001m='soviet-001m.shp'
file500k='soviet-500k.shp'
file200k='soviet-200k.shp'
file100k='soviet-100k.shp'
file050k='soviet-050k.shp'

# output driver
driver = ogr.GetDriverByName('ESRI Shapefile')

driver.DeleteDataSource(file001m)
driver.DeleteDataSource(file500k)
driver.DeleteDataSource(file200k)
driver.DeleteDataSource(file100k)
driver.DeleteDataSource(file050k)

ds001 = driver.CreateDataSource(file001m)
ds500 = driver.CreateDataSource(file500k)
ds200 = driver.CreateDataSource(file200k)
ds100 = driver.CreateDataSource(file100k)
ds050 = driver.CreateDataSource(file050k)

layer001 = ds001.CreateLayer('soviet-001m', geom_type=ogr.wkbPolygon)
layer500 = ds500.CreateLayer('soviet-500k', geom_type=ogr.wkbPolygon)
layer200 = ds200.CreateLayer('soviet-200k', geom_type=ogr.wkbPolygon)
layer100 = ds200.CreateLayer('soviet-100k', geom_type=ogr.wkbPolygon)
layer050 = ds200.CreateLayer('soviet-050k', geom_type=ogr.wkbPolygon)

# create the LABEL field as text
fd = ogr.FieldDefn('LABEL', ogr.OFTString)
fd.SetWidth(12)
fd.SetPrecision(0)

layer001.CreateField(fd)
layer500.CreateField(fd)
layer200.CreateField(fd)
layer100.CreateField(fd)
layer050.CreateField(fd)

# stepping at 1:1.000.000 level
# 6 degree -180 <-> +180
# 4 degree -90 <-> -90
deltax = 6
deltay = 4

iX=-180
while (iX < 180):
    iY=-88
    while (iY <= 88):
        xpos = (iX+180)/deltax + 1
        if iY < 0:
         ypos = (90/deltay) - ((iY+90)/deltay) - 1
         LABEL1M='x%s-%s' % (alphabet[ypos],xpos)
        if iY >= 0:
         ypos = ((90/deltay) - ((iY+90)/deltay)) * -1
         LABEL1M='%s-%s' % (alphabet[ypos],xpos)

        lX001=iX
        rX001=iX+deltax
        uY001=iY
        lY001=iY+deltay

        print '001m: %s [%s %s %s %s]' % (LABEL1M,lX001,uY001,rX001,lY001)

        poly001 = ogr.Geometry(type=ogr.wkbPolygon)
        wkt = 'POLYGON ((%f %f, %f %f, %f %f, %f %f))' \
               % (lX001, uY001, rX001, uY001, rX001, lY001, lX001, lY001)
        poly001 = ogr.CreateGeometryFromWkt(wkt)
        f001 = ogr.Feature(feature_def=layer001.GetLayerDefn())
        f001.SetField(0,'%s' % LABEL1M )
        f001.SetGeometryDirectly(poly001)
        layer001.CreateFeature(f001)
        f001.Destroy()

        SUBLABEL500=1
        for divy in range ( 0, 2 ):
          uY500=lY001 - (divy * (deltay / 2))
          lY500=lY001 - ((divy+1) * deltay / 2)
          for divx in range ( 0, 2 ):
            lX500=lX001 + (divx * deltax / 2)
            rX500=lX001 + ( ( divx+1 ) * (deltax / 2))

            LABEL500K='%s-%s' % (LABEL1M,SUBLABEL500)
            print '--500k: %s [%s %s %s %s]' % (LABEL500K,lX500,uY500,rX500,lY500)

            poly500 = ogr.Geometry(type=ogr.wkbPolygon)
            wkt = 'POLYGON ((%f %f, %f %f, %f %f, %f %f))' \
                   % (lX500, uY500, rX500, uY500, rX500, lY500, lX500, lY500)
            poly500 = ogr.CreateGeometryFromWkt(wkt)
            f500 = ogr.Feature(feature_def=layer500.GetLayerDefn())
            f500.SetField(0,'%s' % LABEL500K )
            f500.SetGeometryDirectly(poly500)
            layer500.CreateFeature(f500)
            f500.Destroy()
            SUBLABEL500=SUBLABEL500+1

        SUBLABEL200=1
        for divy in range ( 0, 6 ):
          uY200=lY001 - (divy * (deltay / 6))
          lY200=lY001 - ((divy+1) * deltay / 6)
          for divx in range ( 0, 6 ):
            lX200=lX001 + (divx * deltax / 6)
            rX200=lX001 + ( ( divx+1 ) * (deltax / 6))

            LABEL200K='%s-%s' % (LABEL1M,SUBLABEL200)
            print '----200k: %s [%s %s %s %s]' % (LABEL200K,lX200,uY200,rX200,lY200)

            poly200 = ogr.Geometry(type=ogr.wkbPolygon)
            wkt = 'POLYGON ((%f %f, %f %f, %f %f, %f %f))' \
                   % (lX200, uY200, rX200, uY200, rX200, lY200, lX200, lY200)
            poly200 = ogr.CreateGeometryFromWkt(wkt)
            f200 = ogr.Feature(feature_def=layer200.GetLayerDefn())
            f200.SetField(0,'%s' % LABEL200K )
            f200.SetGeometryDirectly(poly200)
            layer200.CreateFeature(f200)
            f200.Destroy()
            SUBLABEL200=SUBLABEL200+1

        SUBLABEL100=1
        for divy in range ( 0, 12 ):
          uY100=float(lY001) - (float(divy) * (float(deltay) / 12))
          lY100=float(lY001) - ((float(divy+1)) * float(deltay) / 12)
          for divx in range ( 0, 12 ):
            lX100=float(lX001) + (float(divx) * float(deltax) / 12)
            rX100=float(lX001) + ( ( float(divx+1) ) * (float(deltax) / 12))
            LABEL100K='%s-%s' % (LABEL1M,SUBLABEL100)
            print '------100k: %s [%s %s %s %s]' % (LABEL100K,lX100,uY100,rX100,lY100)

            poly100 = ogr.Geometry(type=ogr.wkbPolygon)
            wkt = 'POLYGON ((%f %f, %f %f, %f %f, %f %f))' \
                   % (lX100, uY100, rX100, uY100, rX100, lY100, lX100, lY100)
            poly100 = ogr.CreateGeometryFromWkt(wkt)
            f100 = ogr.Feature(feature_def=layer100.GetLayerDefn())
            f100.SetField(0,'%s' % LABEL100K )
            f100.SetGeometryDirectly(poly100)
            layer100.CreateFeature(f100)
            f100.Destroy()

            SUBLABEL050=1
            for divy in range ( 0, 2 ):
              uY050=float(lY100) - (float(divy) * (float(deltay) / 24))
              lY050=float(lY100) - ((float(divy+1)) * float(deltay) / 24)
              for divx in range ( 0, 2 ):
                lX050=float(lX100) + (float(divx) * float(deltax) / 24)
                rX050=float(lX100) + ( ( float(divx+1) ) * (float(deltax) / 24))

                LABEL050K='%s-%s' % (LABEL100K,SUBLABEL050)
                print '--------050k: %s [%s %s %s %s]' % (LABEL050K,lX050,uY050,rX050,lY050)

                poly050 = ogr.Geometry(type=ogr.wkbPolygon)
                wkt = 'POLYGON ((%f %f, %f %f, %f %f, %f %f))' \
                       % (lX050, uY050, rX050, uY050, rX050, lY050, lX050, lY050)
                poly050 = ogr.CreateGeometryFromWkt(wkt)
                f050 = ogr.Feature(feature_def=layer050.GetLayerDefn())
                f050.SetField(0,'%s' % LABEL050K )
                f050.SetGeometryDirectly(poly050)
                layer050.CreateFeature(f050)
                f050.Destroy()
                SUBLABEL050=SUBLABEL050+1

            SUBLABEL100=SUBLABEL100+1


        iY=iY+deltay
    iX=iX+deltax