geometry.py 4.47 KB
 Daniel Scheffler committed Aug 07, 2016 1 ``````# -*- coding: utf-8 -*- `````` Daniel Scheffler committed Oct 11, 2016 2 ``````__author__='Daniel Scheffler' `````` Daniel Scheffler committed Apr 20, 2016 3 `````` `````` Daniel Scheffler committed Oct 11, 2016 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 `````` Daniel Scheffler committed Sep 27, 2016 19 `````` `````` Daniel Scheffler committed Oct 11, 2016 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 committed Oct 07, 2016 23 ``````from py_tools_ds.ptds import GeoArray `````` 24 25 `````` `````` Daniel Scheffler committed Apr 20, 2016 26 27 28 29 30 `````` def angle_to_north(XY): """Calculates the angle between the lines [origin:[0,0],north:[0,1]] and `````` Daniel Scheffler committed Oct 11, 2016 31 32 `````` [origin:[0,0],pointXY:[X,Y]] in clockwise direction. Returns values between 0 and 360 degrees. """ `````` Daniel Scheffler committed Apr 20, 2016 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 committed Oct 07, 2016 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 committed Oct 07, 2016 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 committed Oct 07, 2016 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 committed Oct 07, 2016 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 committed Oct 12, 2016 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 committed Oct 07, 2016 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 committed Oct 07, 2016 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 committed Oct 07, 2016 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 committed Oct 07, 2016 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``