# encoding: utf-8
import numpy
import os 
from osgeo import gdal,osr
filePath = "Z:\\ucheba\\BANK_MAP\\Penza\\penza_120.mtw"
dataset = gdal.Open(filePath, gdal.GA_ReadOnly)
print(type(dataset))
print "Driver short name", dataset.GetDriver().ShortName
print("----------------------------------------------------------------")
print "Driver long name", dataset.GetDriver().LongName
print("----------------------------------------------------------------")
print "Raster size", dataset.RasterXSize, "x", dataset.RasterYSize
print("----------------------------------------------------------------")
print "Number of bands", dataset.RasterCount
print("----------------------------------------------------------------")
old_coordinate_system = osr.SpatialReference()
old_coordinate_system.ImportFromWkt(dataset.GetProjectionRef())
print "Projection", old_coordinate_system
print("----------------------------------------------------------------")
print "Geo transform", dataset.GetGeoTransform()
print("----------------------------------------------------------------")
width = dataset.RasterXSize
height = dataset.RasterYSize
gt = dataset.GetGeoTransform()
minx = gt[0]
miny = gt[3] + width*gt[4] + height*gt[5] 
maxx = gt[0] + width*gt[1] + height*gt[2]
maxy = gt[3] 
print "X min",minx
print "Y min",miny
print "X max",maxx
print "Y max",maxy
print("----------------------------------------------------------------")

# create the new coordinate system
wgs84_wkt = """
GEOGCS["Pulkovo 1942",
        DATUM["Pulkovo_1942",
            SPHEROID["Krassowsky 1940",6378245,298.3,
                AUTHORITY["EPSG","7024"]],
            TOWGS84[23.92,-141.27,-80.9,0,0.35,0.82,-0.12],
            AUTHORITY["EPSG","6284"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4284"]]"""
new_cs = osr.SpatialReference()
new_cs .ImportFromWkt(wgs84_wkt)

# create a transform object to convert between coordinate systems
transform = osr.CoordinateTransformation(old_coordinate_system,new_cs)
print "transform",transform

#get the point to transform, pixel (0,0) in this case
width = dataset.RasterXSize
print "width",width
height = dataset.RasterYSize
print "height",height
gt = dataset.GetGeoTransform()
min_X = gt[0]
min_Y = gt[3] + width*gt[4] + height*gt[5] 
print(min_X)
print(min_Y)
#get the coordinates in lat long
latlong = transform.TransformPoint(min_X,min_Y)
latlong_1 = transform.TransformPoint(maxx,maxy)
print u"Нижний левый угол квардрата",latlong 
print u"Правый верхний угол квардрата",latlong_1
print("----------------------------------------------------------------")
raster = dataset.ReadAsArray()
#print(raster)
#print("----------------------------------------------------------------")
gdalBand = dataset.GetRasterBand( 1 )
#band = gdalBand.ReadAsArray()
#print(band)
#print("----------------------------------------------------------------")
line = gdalBand.ReadAsArray( 0, 0, 1781, 1994 )#1781, 1994
print(line)
print("----------------------------------------------------------------")

#places_list = line
#with open('W:\\Desktop\\TXT.txt', 'w') as filehandle:  
#    filehandle.writelines("%s\n" % place for place in places_list)
#ds = gdal.Translate('output.txt', ds)
print("Hello World")
input('Press ENTER to exit')