#!/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 MGI - Militar Geographische Institut (Former Yugoslavia & Austria)
# grid at 100k/050k/025k in booth in standard lonlat degree coordinates and
# steppings as generic EPSG:4312 describes.
#
# 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

# stepping at 1:100k level
delta = .5

# output files
file100k='mgi-100k.shp'
file050k='mgi-050k.shp'
file025k='mgi-025k.shp'

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

driver.DeleteDataSource(file100k)
driver.DeleteDataSource(file050k)
driver.DeleteDataSource(file025k)

ds100 = driver.CreateDataSource(file100k)
ds050 = driver.CreateDataSource(file050k)
ds025 = driver.CreateDataSource(file025k)

layer100 = ds100.CreateLayer('mgi-100k', geom_type=ogr.wkbPolygon)
layer050 = ds050.CreateLayer('mgi-050k', geom_type=ogr.wkbPolygon)
layer025 = ds025.CreateLayer('mgi-025k', geom_type=ogr.wkbPolygon)

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

layer100.CreateField(fd)
layer050.CreateField(fd)
layer025.CreateField(fd)

# IME (.yu) = NUME (.ro)
# IME (.yu) = LABEL (.en)

# generate MGI coverage of 100k
IME100=1
for iY100 in range ( 1, 51 ):
    uY100 = 49.5 - ( iY100 * delta )
    lY100 = 49.5 - ( ( iY100 +1 ) * delta )
    for iX100 in range ( 1, 51 ):
        lX100 = 5.5 + ( iX100 * delta )
        rX100 = 5.5 + ( ( iX100+1 ) * delta )

        print '100k: %s: [%s %s] [%s %s]' % (IME100,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' % IME100 )
        f100.SetGeometryDirectly(poly100)
        layer100.CreateFeature(f100)
        f100.Destroy()

        # lower/generate MGI coverage to 50k
        IME50=1
        for iY50 in range ( 0 , 2 ):
            uY50=uY100 - (iY50 * (delta / 2))
            lY50=uY100 - ((iY50+1) * delta / 2)
            for iX50 in range ( 0 , 2 ):
                lX50=lX100 + (iX50 * delta / 2)
                rX50=lX100 + ( ( iX50+1 ) * (delta / 2))

                print '--050k: %s-%s: [%s %s] [%s %s]' % (IME100,IME50,lX50,uY50,rX50,lY50)
                poly050 = ogr.Geometry(type=ogr.wkbPolygon)
                wkt = 'POLYGON ((%f %f, %f %f, %f %f, %f %f))' \
                       % (lX50, uY50, rX50, uY50, rX50, lY50, lX50, lY50)
                poly050 = ogr.CreateGeometryFromWkt(wkt)
                f050 = ogr.Feature(feature_def=layer050.GetLayerDefn())
                f050.SetField(0,'%s-%s' % (IME100,IME50) )
                f050.SetGeometryDirectly(poly050)
                layer050.CreateFeature(f050)
                f050.Destroy()


                # lower/generate MGI coverage to 25k
                IME25=1
                for iY25 in range ( 0 , 2 ):
                    uY25=uY50 - (iY25 * (delta / 4))
                    lY25=uY50 - ((iY25+1) * delta / 4)
                    for iX25 in range ( 0 , 2 ):
                        lX25=lX50 + (iX25 * delta / 4)
                        rX25=lX50 + ( ( iX25+1 ) * (delta / 4))

                        print '----025k: %s-%s-%s: [%s %s] [%s %s]' % (IME100,IME50,IME25,lX25,uY25,rX25,lY25)
                        poly025 = ogr.Geometry(type=ogr.wkbPolygon)
                        wkt = 'POLYGON ((%f %f, %f %f, %f %f, %f %f))' \
                               % (lX25, uY25, rX25, uY25, rX25, lY25, lX25, lY25)
                        poly025 = ogr.CreateGeometryFromWkt(wkt)
                        f025 = ogr.Feature(feature_def=layer025.GetLayerDefn())
                        f025.SetField(0,'%s-%s-%s' % (IME100,IME50,IME25) )
                        f025.SetGeometryDirectly(poly025)
                        layer025.CreateFeature(f025)
                        f025.Destroy()

                        IME25=IME25+1

                IME50=IME50+1

        IME100=IME100+1