geometry.py 4.47 KB
Newer Older
1
# -*- coding: utf-8 -*-
2
__author__='Daniel Scheffler'
3

4
import warnings
5
6
7
import sys

# custom
8
import numpy  as np
9
from geopandas import GeoDataFrame
10
11
12
13
14
15
16
17
18

try:
    import gdal
    import osr
    import ogr
except ImportError:
    from osgeo import gdal
    from osgeo import osr
    from osgeo import ogr
19

20
21
22
# internal modules
from py_tools_ds.ptds.geo.coord_calc       import calc_FullDataset_corner_positions
from py_tools_ds.ptds.geo.coord_trafo      import pixelToMapYX, imYX2mapYX
Daniel Scheffler's avatar
GEO:    
Daniel Scheffler committed
23
from py_tools_ds.ptds                      import GeoArray
24
25


26
27
28
29
30



def angle_to_north(XY):
    """Calculates the angle between the lines [origin:[0,0],north:[0,1]] and
31
32
     [origin:[0,0],pointXY:[X,Y]] in clockwise direction. Returns values between 0 and 360 degrees.
     """
33
34
35
36
37
    XY    = np.array(XY)
    XYarr = XY if len(XY.shape)==2 else XY.reshape((1,2))
    return np.abs(np.degrees(np.arctan2(XYarr[:,1],XYarr[:,0])-np.pi/2)%360)


38
def get_true_corner_mapXY(fPath_or_geoarray, bandNr=1, noDataVal=None, mp=1, v=0, q=0):
Daniel Scheffler's avatar
GEO:    
Daniel Scheffler committed
39
40
41
42
    geoArr    = GeoArray(fPath_or_geoarray) if not isinstance(fPath_or_geoarray,GeoArray) else fPath_or_geoarray

    rows,cols = geoArr.shape[:2]
    gt, prj   = geoArr.geotransform, geoArr.projection
43

Daniel Scheffler's avatar
GEO:    
Daniel Scheffler committed
44
    assert gt and prj, 'GeoTransform an projection must be given for calculation of LonLat corner coordinates.'
45
46
47
48
49

    mask_1bit = np.zeros((rows,cols),dtype='uint8') # zeros -> image area later overwritten by ones

    if noDataVal is None:
        mask_1bit[:,:] = 1
50
    elif noDataVal=='ambiguous':
51
        warnings.warn("No data value could not be automatically detected. Thus the matching window used for shift "
52
              "calculation had to be centered in the middle of the overlap area without respecting no data values. "
53
54
55
              "To avoid this provide the correct no data values for reference and shift image via '-nodata'")
        mask_1bit[:,:] = 1
    else:
Daniel Scheffler's avatar
GEO:    
Daniel Scheffler committed
56
        band_data = geoArr[bandNr-1] # TODO implement gdal_ReadAsArray_mp (reading in multiprocessing)
57
58
59
60
        mask_1bit[band_data!=noDataVal] = 1

    if v: print('detected no data value',noDataVal)

Daniel Scheffler's avatar
GEO:    
Daniel Scheffler committed
61
62
    try:
        corner_coords_YX = calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=False, algorithm='shapely')
63
64
65
66
    except Exception:
        if v:
            warnings.warn("\nCalculation of corner coordinates failed within algorithm 'shapely' (Exception: %s)."
                          " Using algorithm 'numpy' instead." %sys.exc_info()[1])
67
        corner_coords_YX = calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=False, algorithm='numpy') # FIXME numpy algorithm returns wrong values for S2A_OPER_MSI_L1C_TL_SGS__20160608T153121_A005024_T33UUU_B03.jp2 (Hannes)
68

Daniel Scheffler's avatar
CoReg:    
Daniel Scheffler committed
69
70
71
    if len(corner_coords_YX)==4: # this avoids shapely self intersection
        corner_coords_YX = list(np.array(corner_coords_YX)[[0,1,3,2]]) # UL, UR, LL, LR => UL, UR, LR, LL

72
73
74
75
76
77
78
79
    # check if enough unique coordinates have been found
    if not len(GeoDataFrame(corner_coords_YX).drop_duplicates().values)>=3:
        if not q:
            warnings.warn('\nThe algorithm for automatically detecting the actual image coordinates did not find '
                          'enough unique corners. Using outer image corner coordinates instead.')
        corner_coords_YX = ((0, 0), (0, cols-1), (rows-1, 0), (rows-1, cols-1))

    # check if all points are unique
Daniel Scheffler's avatar
GEO:    
Daniel Scheffler committed
80
81
    #all_coords_are_unique = len([UL, UR, LL, LR]) == len(GeoDataFrame([UL, UR, LL, LR]).drop_duplicates().values)
    #UL, UR, LL, LR = (UL, UR, LL, LR) if all_coords_are_unique else ((0, 0), (0, cols-1), (rows-1, 0), (rows-1, cols-1))
82

83
    get_mapYX     = lambda YX: pixelToMapYX(list(reversed(YX)),geotransform=gt,projection=prj)[0]
Daniel Scheffler's avatar
GEO:    
Daniel Scheffler committed
84
    corner_pos_XY = [list(reversed(i)) for i in [get_mapYX(YX) for YX in corner_coords_YX]]
85
86
87
88
89
90
91
92
93
94
95
96
97
    return corner_pos_XY


def get_subset_GeoTransform(gt_fullArr,subset_box_imYX):
    gt_subset = list(gt_fullArr[:]) # copy
    gt_subset[3],gt_subset[0] = imYX2mapYX(subset_box_imYX[0],gt_fullArr)
    return gt_subset


def get_gdalReadInputs_from_boxImYX(boxImYX):
    """Returns row_start,col_start,rows_count,cols_count and assumes boxImYX as [UL_YX,UR_YX,LR_YX,LL_YX)"""
    rS, cS = boxImYX[0]
    clip_sz_x = abs(boxImYX[1][1]-boxImYX[0][1]) # URx-ULx
Daniel Scheffler's avatar
GEO:    
Daniel Scheffler committed
98
    clip_sz_y = abs(boxImYX[0][0]-boxImYX[3][0]) # ULy-LLy
99
100
101
    return cS, rS, clip_sz_x,clip_sz_y


Daniel Scheffler's avatar
GEO:    
Daniel Scheffler committed
102
103
104
105
def get_GeoArrayPosition_from_boxImYX(boxImYX):
    """Returns row_start,row_end,col_start,col_end and assumes boxImYX as [UL_YX,UR_YX,LR_YX,LL_YX)"""
    rS, cS = boxImYX[0] # UL
    rE, cE = boxImYX[2] # LR
106
    return rS, rE, cS, cE