map_info.py 3.4 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"


try:
    from osgeo import osr
except ImportError:
    import osr

from .coord_trafo import transform_utm_to_wgs84


def geotransform2mapinfo(gt,prj):
    """Builds an ENVI map info from given GDAL GeoTransform and Projection (compatible with UTM and LonLat projections).
    :param gt:  GDAL GeoTransform, e.g. (249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0)
    :param prj: GDAL Projection - WKT Format
    :returns:   ENVI map info, e.g. [ UTM , 1 , 1 , 256785.0 , 4572015.0 , 30.0 , 30.0 , 43 , North , WGS-84 ]
    :rtype:     list
    """

    if gt[2]!=0 or gt[4]!=0: # TODO
        raise NotImplementedError('Currently rotated datasets are not supported.')
    srs         = osr.SpatialReference()
    srs.ImportFromWkt(prj)
    Proj4       = [i[1:] for i in srs.ExportToProj4().split()]
    Proj4_proj  = [v.split('=')[1] for i,v in enumerate(Proj4) if '=' in v and v.split('=')[0]=='proj'][0]
    Proj4_ellps = [v.split('=')[1] for i,v in enumerate(Proj4) if '=' in v and v.split('=')[0] in ['ellps','datum']][0]
    proj        = 'Geographic Lat/Lon' if Proj4_proj=='longlat' else 'UTM' if Proj4_proj=='utm' else Proj4_proj
    ellps       = 'WGS-84' if Proj4_ellps=='WGS84' else Proj4_ellps
    UL_X, UL_Y, GSD_X, GSD_Y = gt[0], gt[3], gt[1], gt[5]
    utm2wgs84          = lambda utmX,utmY: transform_utm_to_wgs84(utmX,utmY,srs.GetUTMZone())
    is_UTM_North_South = lambda LonLat: 'North' if LonLat[1] >= 0. else 'South'
    map_info = [proj,1,1,UL_X,UL_Y,GSD_X,abs(GSD_Y),ellps] if proj!='UTM' else \
        ['UTM',1,1,UL_X,UL_Y,GSD_X,abs(GSD_Y),srs.GetUTMZone(),is_UTM_North_South(utm2wgs84(UL_X,UL_Y)),ellps]
    return map_info


def mapinfo2geotransform(map_info):
    # type: (list) -> list
    """Builds GDAL GeoTransform tuple from an ENVI map info.
    :param map_info: ENVI map info (list), e.g. ['UTM', 1, 1, 192585.0, 5379315.0, 30.0, 30.0, 41, 'North', 'WGS-84']
    :returns:        GDAL GeoTransform, e.g. [249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0]
    :rtype:          list
    """
    # FIXME rotated datasets are currently not supported -> function must return rotation at gt[2] and gt[4]
    if map_info[1]==1 and map_info[2]==1:
        ULmapX, ULmapY = float(map_info[3]),float(map_info[4])
    else:
        ULmapX = float(map_info[3])-(float(map_info[1])*float(map_info[5])-float(map_info[5]))
        ULmapY = float(map_info[4])+(float(map_info[2])*float(map_info[6])-float(map_info[6]))
    assert isinstance(map_info,list) and len(map_info)==10,\
        "The map info argument has to be a list of length 10. Got %s." %map_info
    return [ULmapX,float(map_info[5]),0.,ULmapY,0.,-float(map_info[6])]

def get_corner_coordinates(gdal_ds=None, gt=None, cols=None, rows=None):
    """Returns (ULxy, LLxy, LRxy, URxy) in the same coordinate units like the given geotransform."""
    assert gdal_ds or (gt and cols and rows), \
        "GEOP.get_corner_coordinates: Missing argument! Please provide either 'gdal_ds' or 'gt', 'cols' AND 'rows'."
    gdal_ds_GT = gdal_ds.GetGeoTransform() if gdal_ds else gt
    ext=[]
    xarr=[0,gdal_ds.RasterXSize if gdal_ds else cols]
    yarr=[0,gdal_ds.RasterYSize if gdal_ds else rows]
    for px in xarr:
        for py in yarr:
            x=gdal_ds_GT[0]+(px*gdal_ds_GT[1])+(py*gdal_ds_GT[2])
            y=gdal_ds_GT[3]+(px*gdal_ds_GT[4])+(py*gdal_ds_GT[5])
            ext.append([x,y])
        yarr.reverse()
    gdal_ds_GT = None
    return ext