Перевод координат

[Starik]

Здравствуйте, уважаемые.
У меня проблема. Хочу перевести координаты из географических в координаты пиксела на растре через конвертацию в UTM.
Предположим есть 2 привязаных растра в GeoTIFF. Gdal выдает на них вот-такую информацию:
Растр №1:

Код: Выделить всё

Projection GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]

Geo transform (5363415.843844333, 7.007168368629773, 0.0, 7223749.626219502, 0.0, -7.007168368629773)

Растр №2:

Код: Выделить всё

Projection PROJCS["unnamed",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223560493,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","3395"]]

Geo transform (5009177.638276306, 108.2969305156259, 0.0, 7522996.349775525, 0.0, -108.29693051562595)

У меня получилась вот такая функция для этого перевода:

Код: Выделить всё

from osgeo import gdal, osr

def getXY(gdalData, E, N):
    """
    Делает из координаты типа 49, 54.9 - координату в пикселах на изображении
    """
    # get the existing coordinate system
    new_cs= osr.SpatialReference()
    new_cs.ImportFromWkt(gdalData.GetProjectionRef())

    # create the new coordinate system
    wgs84_wkt = """
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]]"""
    old_cs = osr.SpatialReference()
    old_cs.ImportFromWkt(wgs84_wkt)

    # create a transform object to convert between coordinate systems
    transform = osr.CoordinateTransformation(old_cs,new_cs)
    gt = gdalData.GetGeoTransform()

    width, height = gdalData.RasterXSize, gdalData.RasterYSize
    x, y, z = transform.TransformPoint(E, N)

    #print x, y, z, E, N,

    x = int((x - gt[0]) / gt[1])
    y = int((y - gt[3]) / gt[5])

    if x < 0 or y < 0 or x > width or y > height:
        return False

    return x, y

Код: Выделить всё

gdalData = gdal.Open(inputF)

, где inputF - путь к файлу
E = Долгота
N = Широта

На растре №2 работает идеально, а на растре №1 возвращает те же E и N.
Я понимаю что где то в определении координатных систем для трансформации косяк, но опыта у меня маловато. Мне хоть куда копать узнать бы.
Заранее благодарен.

[Starik]

Эх, специалисты, где же вы где. Задачу решил, но так и не понял как это все работает, но главное что работает.

Код: Выделить всё

def getXY(gdalData, E, N):

    CS = [

        'PROJCS["unnamed",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223560493,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","3395"]]',

        'GEOGCS["WGS 84",DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.01745329251994328, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4326"]]',

        gdalData.GetProjectionRef()

    ]

    inpCs, outCs = osr.SpatialReference(), osr.SpatialReference()
    gt = gdalData.GetGeoTransform()
    width, height = gdalData.RasterXSize, gdalData.RasterYSize

    for i in CS:
        for j in CS:
            outCs.ImportFromWkt(i)
            inpCs.ImportFromWkt(j)
            transform = osr.CoordinateTransformation(inpCs, outCs)
            x, y, z = transform.TransformPoint(E, N)

            x = int((x - gt[0]) / gt[1])
            y = int((y - gt[3]) / gt[5])

            if width >= x >= 0 and height >= y >= 0:
                return x, y

    return False