1# -*- coding: utf-8 -*- 

2 

3# gms_preprocessing, spatial and spectral homogenization of satellite remote sensing data 

4# 

5# Copyright (C) 2020 Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de) 

6# 

7# This software was developed within the context of the GeoMultiSens project funded 

8# by the German Federal Ministry of Education and Research 

9# (project grant code: 01 IS 14 010 A-C). 

10# 

11# This program is free software: you can redistribute it and/or modify it under 

12# the terms of the GNU General Public License as published by the Free Software 

13# Foundation, either version 3 of the License, or (at your option) any later version. 

14# Please note the following exception: `gms_preprocessing` depends on tqdm, which 

15# is distributed under the Mozilla Public Licence (MPL) v2.0 except for the files 

16# "tqdm/_tqdm.py", "setup.py", "README.rst", "MANIFEST.in" and ".gitignore". 

17# Details can be found here: https://github.com/tqdm/tqdm/blob/master/LICENCE. 

18# 

19# This program is distributed in the hope that it will be useful, but WITHOUT 

20# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 

21# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more 

22# details. 

23# 

24# You should have received a copy of the GNU Lesser General Public License along 

25# with this program. If not, see <http://www.gnu.org/licenses/>. 

26 

27import shapely 

28from py_tools_ds.geo.coord_trafo import reproject_shapelyGeometry 

29from shapely.geometry import Polygon 

30 

31from gms_preprocessing.options.config import GMS_config as CFG 

32from gms_preprocessing.misc import database_tools as DB_T 

33from gms_preprocessing.misc import helper_functions as HLP_F 

34from gms_preprocessing.algorithms import geoprocessing as GEOP 

35 

36__author__ = 'Daniel Scheffler' 

37 

38 

39class MGRS_tile(object): 

40 def __init__(self, tile_ID=''): 

41 """ 

42 

43 :param tile_ID: <str> 5 digit tile ID, e.g. 32UUU 

44 """ 

45 self._tile_ID = '' 

46 self._geom_wkb = '' 

47 self._poly_lonlat = Polygon() 

48 

49 if tile_ID: 

50 self.tile_ID = tile_ID 

51 

52 @property 

53 def tile_ID(self): 

54 return self._tile_ID 

55 

56 @tile_ID.setter 

57 def tile_ID(self, tile_ID): 

58 assert isinstance(tile_ID, str) and len(tile_ID) == 5, \ 

59 "'tile_ID' must be a 5 digit string code. Got %s." % tile_ID 

60 self._tile_ID = tile_ID 

61 

62 res = DB_T.get_info_from_postgreSQLdb(CFG.conn_database, 'mgrs_tiles', ['geom'], 

63 {'grid1mil': self.grid1mil, 'grid100k': self.grid100k}, timeout=20000) 

64 assert res, "The tile ID '%s' does not exist in the database." % tile_ID 

65 

66 self.geom_wkb = res[0][0] 

67 

68 @property 

69 def grid1mil(self): 

70 return self.tile_ID[:3] 

71 

72 @property 

73 def grid100k(self): 

74 return self.tile_ID[-2:] 

75 

76 @property 

77 def UTMzone(self): 

78 return int(self.tile_ID[:2]) 

79 

80 @property 

81 def EPSG(self): 

82 is_south = self.poly_lonlat.centroid.xy[1][0] < 0 

83 return int(('327' if is_south else '326') + str(self.UTMzone)) 

84 

85 @property 

86 def geom_wkb(self): 

87 return self._geom_wkb 

88 

89 @geom_wkb.setter 

90 def geom_wkb(self, geom_wkb): 

91 # FIXME should check whether geom_wkb is a valid geometry according to db 

92 self._geom_wkb = geom_wkb 

93 

94 @property 

95 def poly_lonlat(self): 

96 return shapely.wkb.loads(self.geom_wkb, hex=True) 

97 

98 @property 

99 def poly_utm(self): 

100 return reproject_shapelyGeometry(self.poly_lonlat, 4326, self.EPSG) 

101 

102 def poly_specPrj(self, prj): 

103 """Returns a shapely.Polygon in a specific projection. 

104 

105 :param prj: <str> WKT string of the target projection 

106 """ 

107 return reproject_shapelyGeometry(self.poly_lonlat, 4326, prj) 

108 

109 def get_bounds(self, prj=None): 

110 return self.poly_lonlat.bounds if not prj else self.poly_specPrj(prj).bounds # xmin, ymin, xmax, ymax 

111 

112 def to_image_bounds(self, im_prj, im_gt, arr_shape, pixbuffer=0, ensure_valid_coords=True): 

113 xmin, xmax, ymin, ymax = HLP_F.get_arrSubsetBounds_from_shapelyPolyLonLat(arr_shape, self.poly_lonlat, im_gt, 

114 im_prj, pixbuffer, 

115 ensure_valid_coords) 

116 return xmin, xmax, ymin, ymax 

117 

118 def clip_array_using_mgrsBounds(self, array, im_prj, im_gt, nodataVal=0, pixbuffer=0): 

119 """ 

120 

121 :param array: 

122 :param im_prj: 

123 :param im_gt: 

124 :param nodataVal: 

125 :param pixbuffer: <float> an optional buffer size (image pixel units) 

126 """ 

127 

128 buffMap = im_gt[1] * pixbuffer 

129 mgrs_bounds = self.poly_specPrj(im_prj).buffer(buffMap).bounds 

130 tgt_arr, tgt_gt, im_prj = GEOP.clip_array_using_mapBounds(array, mgrs_bounds, im_prj, im_gt, nodataVal) 

131 

132 return tgt_arr, tgt_gt, im_prj