Commit d9aeb8bd authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Minor code improvements.

parent 8400a748
Pipeline #1269 passed with stage
in 16 minutes and 23 seconds
......@@ -237,7 +237,8 @@ class AtmCorr(object):
# [setattr(L1C_obj, 'AtmCorr', self) for L1C_obj in self.inObjs] # too big for serialization
if not re.search('Sentinel-2', self.inObjs[0].satellite, re.I):
warnings.warn('Calculation of acquisition geometry arrays is currently only validated for Sentinel-2!')
self.logger.warning('Calculation of acquisition geometry arrays is currently only validated for '
'Sentinel-2!')
# validation possible by comparing S2 angles provided by ESA with own angles
@property
......
......@@ -48,8 +48,6 @@ __author__ = 'Daniel Scheffler', 'Robert Behling'
class GEOPROCESSING(object):
"""****CREATE OBJECT ****************************************************"""
"""****OBJECT CONSTRUCTOR**************************************************"""
def __init__(self, geodata, logger, workspace=None, subset=None, v=None):
self.logger = logger
......@@ -210,8 +208,8 @@ class GEOPROCESSING(object):
self.Prop = {'DataProp': self.DataProp, 'DriverProp': self.DriverProp, 'GeoProp': self.GeoProp}
def gdalinfo(self):
"""----METHOD_1----------------------------------------------------------
get properties of the Inputdatasets via gdalinfo
"""get properties of the Inputdatasets via gdalinfo
(die Infos die er hier ausschreibt in subrocess.Popen können einer variable als Text übergeben werden. Aus
diesem Text kann ich dann die Infos als Attribute ausschreiben.
als ersatz für die jetzige Attributerzeugung
......@@ -223,6 +221,7 @@ class GEOPROCESSING(object):
def add_GeoTransform_Projection_using_MetaData(self, CornerTieP_LonLat, CS_EPSG=None, CS_TYPE=None, CS_DATUM=None,
CS_UTM_ZONE=None, gResolution=None, shape_fullArr=None):
""" Method assumes that north is up. Map rotations are not respected.
:param CornerTieP_LonLat:
:param CS_EPSG:
:param CS_TYPE:
......@@ -282,6 +281,7 @@ class GEOPROCESSING(object):
dst_CS='UTM', dst_CS_datum='WGS84', mode='GDAL', use_workspace=True,
inFill=None):
"""Warp image to new Projection.
:param use_inherent_GCPs:
:param TieP: Corner Tie Points - always LonLat.
:param dst_EPSG_code: EPSG-Code defines LonLat or UTM coordinates.
......@@ -442,11 +442,14 @@ class GEOPROCESSING(object):
def get_corner_coordinates(self, targetProj):
"""Returns corner coordinates of the entire GEOP object in lon/lat or UTM.
ATTENTION: coordinates represent PIXEL CORNERS:
UL=UL-coordinate of (0,0)
UR=UR-coordinate of (0,self.cols-1)
=> lonlat_arr always contains UL-coordinate of each pixel
:param targetProj: 'LonLat' or 'UTM' """
ATTENTION: coordinates represent PIXEL CORNERS:
UL=UL-coordinate of (0,0)
UR=UR-coordinate of (0,self.cols-1)
=> lonlat_arr always contains UL-coordinate of each pixel
:param targetProj: 'LonLat' or 'UTM'
"""
UL_LonLat = self.pixelToWorldCoord_using_geotransform_and_projection((0, 0), targetProj)
# self.cols statt self.cols-1, um Koordinate am rechten Pixelrand zu berechnen:
UR_LonLat = self.pixelToWorldCoord_using_geotransform_and_projection((0, self.cols), targetProj)
......@@ -460,10 +463,10 @@ class GEOPROCESSING(object):
def calc_mask_data_nodata(self, array=None, custom_nodataVal=-9999):
"""Berechnet den Bildbereich, der von allen Kanälen abgedeckt wird.
:param array: input numpy array to be used for mask calculation (otherwise read from disk)
:param custom_nodataVal:
"""
nodataVal = get_outFillZeroSaturated(np.int16)[0] if custom_nodataVal is None else custom_nodataVal
in_arr = array if array is not None else np.swapaxes(np.swapaxes(gdalnumeric.LoadFile(self.desc), 0, 2), 0, 1)
return np.all(np.where(in_arr == nodataVal, 0, 1), axis=2)
......@@ -485,8 +488,8 @@ class GEOPROCESSING(object):
return transform_wgs84_to_utm(XWorld, YWorld) # [HW,RW]
def tondarray(self, direction=1, startpix=None, extent=None, UL=None, LR=None, v=0):
"""----METHOD_2----------------------------------------------------------
convert gdalobject to 3dimensional ndarray stack ([x,y,z])
"""Convert gdalobject to 3dimensional ndarray stack ([x,y,z]).
:param direction: 1: shape: [bands, rows, cols]
-> same as theres read envi structure ([band1],[band2],..[bandn])
2: shape: [rows, bands, cols]
......@@ -499,11 +502,6 @@ class GEOPROCESSING(object):
:param LR: [x,y] map coordinates -> LR map coordinates for subset generation
:param v:
"""
if v:
print("\n--------GEOPROCESSING--------\n##Method##"
"\n**tondarray**")
# -- test of given inputparameter
test = [startpix is not None, extent is not None, UL is not None, LR is not None]
......@@ -726,9 +724,10 @@ class GEOPROCESSING(object):
def ndarray2gdal(ndarray, outPath=None, importFile=None, direction=1, GDAL_Type="ENVI", geotransform=None,
projection=None, v=0):
"""----FUNCTION_6----------------------------------------------------------
converts a numpy array to a Georasterdataset (default envi bsq)
"""Converts a numpy array to a Georasterdataset (default envi bsq).
bis jetzt nur georeferenzierung und projectionszuweisung durch import aus anderem file möglich
:param ndarray: [bands, rows, cols]
:param outPath: Path were to store the outputFile
:param importFile: path of a file, where to import the projection and mapinfos
......@@ -808,14 +807,12 @@ def ndarray2gdal(ndarray, outPath=None, importFile=None, direction=1, GDAL_Type=
def convertGdalNumpyDataType(dType):
"""----FUNCTION_7----------------------------------------------------------
convertGdalNumpyDataType
"""convertGdalNumpyDataType
input:
dType: GDALdataType string or numpy dataType
output:
corresponding dataType
"""
# dictionary to translate GDAL data types (strings) in corresponding numpy data types
dTypeDic = {"Byte": np.uint8, "UInt16": np.uint16, "Int16": np.int16, "UInt32": np.uint32, "Int32": np.int32,
"Float32": np.float32, "Float64": np.float64, "GDT_UInt32": np.uint32}
......@@ -860,7 +857,6 @@ def get_common_extent(list_extentcoords, alg='outer', return_box=True):
:param return_box: whether to return a coordinate box/envelope of the output Polygon
:return: [ULxy,URxy,LRxy,LLxy]
"""
if alg != 'outer':
raise NotImplementedError # TODO apply shapely intersection
......@@ -905,8 +901,10 @@ def adjust_acquisArrProv_to_shapeFullArr(arrProv, shapeFullArr, meshwidth=1, sub
def DN2Rad(ndarray, offsets, gains, inFill=None, inZero=None, inSaturated=None, cutNeg=True):
# type: (np.ndarray,list,list,int,int,int,bool) -> np.ndarray
"""Convert DN to Radiance [W * m-2 * sr-1 * micrometer-1]
!!! InputGains and Offsets should be in [W * m-2 * sr-1 * micrometer-1]
"""Convert DN to Radiance [W * m-2 * sr-1 * micrometer-1].
NOTE: InputGains and Offsets should be in [W * m-2 * sr-1 * micrometer-1]!!!
:param ndarray: <np.ndarray> array of DNs to be converted into radiance
:param offsets: [W * m-2 * sr-1 * micrometer-1]:
list that includes the offsets of the individual rasterbands [offset_band1, offset_band2,
......@@ -919,7 +917,6 @@ def DN2Rad(ndarray, offsets, gains, inFill=None, inZero=None, inSaturated=None,
:param inSaturated: pixelvalues allocated to saturated pixels
:param cutNeg: cutNegvalues -> all negative values set to 0
"""
assert isinstance(offsets, list) and isinstance(gains, list), \
"Offset and Gain parameters have to be provided as two lists containing gains and offsets for \
each band in ascending order. Got offsets as type '%s' and gains as type '%s'." % (type(offsets), type(gains))
......@@ -953,6 +950,7 @@ def DN2TOARef(ndarray, offsets, gains, irradiances, zenith, earthSunDist,
inFill=None, inZero=None, inSaturated=None, cutNeg=True):
# type: (np.ndarray,list,list,list,float,float,int,int,int,bool) -> np.ndarray
"""Converts DN data to TOA Reflectance.
:param ndarray: <np.ndarray> array of DNs to be converted into TOA reflectance
:param offsets: list: of offsets of each rasterband [W * m-2 * sr-1 * micrometer-1]
[offset_band1, offset_band2, ... ,offset_bandn] or optional as number if Dataset has
......@@ -970,7 +968,6 @@ def DN2TOARef(ndarray, offsets, gains, irradiances, zenith, earthSunDist,
:param cutNeg: bool: if true. all negative values turned to zero. default: True
:return: Int16 TOA_Reflectance in [0-10000]
"""
assert isinstance(offsets, list) and isinstance(gains, list) and isinstance(irradiances, list), \
"Offset, Gain, Irradiance parameters have to be provided as three lists containing gains, offsets and " \
"irradiance for each band in ascending order. Got offsets as type '%s', gains as type '%s' and irradiance as " \
......@@ -1006,6 +1003,7 @@ def TOARad2Kelvin_fastforward(ndarray, K1, K2, emissivity=0.95, inFill=None, inZ
# type: (np.ndarray,list,list,float,int,int,int) -> np.ndarray
"""Convert top-of-atmosphere radiances of thermal bands to temperatures in Kelvin
by applying the inverse of the Planck function.
:param ndarray: <np.ndarray> array of TOA radiance values to be converted into Kelvin
:param K1:
:param K2:
......@@ -1014,7 +1012,6 @@ def TOARad2Kelvin_fastforward(ndarray, K1, K2, emissivity=0.95, inFill=None, inZ
:param inZero:
:param inSaturated:
"""
bands = 1 if len(ndarray.shape) == 2 else ndarray.shape[2]
for arg, argname in zip([K1, K2], ['K1', 'K2']):
assert isinstance(arg[0], float) or isinstance(arg[0], int), "TOARad2Kelvin_fastforward: Expected float or " \
......@@ -1059,7 +1056,6 @@ def DN2DegreesCelsius_fastforward(ndarray, offsets, gains, K1, K2, emissivity=0.
:param inZero:
:param inSaturated:
"""
bands = 1 if len(ndarray.shape) == 2 else ndarray.shape[2]
for arg, argname in zip([offsets, gains, K1, K2], ['Offset', 'Gain', 'K1', 'K2']):
assert isinstance(offsets, list) and isinstance(gains, list), \
......@@ -1121,8 +1117,8 @@ def get_lonlat_coord_array(shape_fullArr, arr_pos, geotransform, projection, mes
:param meshwidth: <int> defines the density of the mesh used for generating the output
(1: full resolution; 10: one point each 10 pixels)
:param nodata_mask: <numpy array>, used for declaring nodata values in the output lonlat array
:param outFill: the value that is assigned to nodata area in the output lonlat array"""
:param outFill: the value that is assigned to nodata area in the output lonlat array
"""
rows, cols, bands, rowStart, rowEnd, colStart, colEnd = get_subsetProps_from_shapeFullArr_arrPos(shape_fullArr,
arr_pos)
......@@ -1170,7 +1166,6 @@ def calc_VAA_using_fullSceneCornerLonLat(fullSceneCornerLonLat, orbit_params):
:param fullSceneCornerLonLat:
:param orbit_params: list of [altitude, inclination, period] => inclination is used as fallback
"""
UL_LonLat, UR_LonLat, LL_LonLat, LR_LonLat = fullSceneCornerLonLat
forward_az_left = pyproj.Geod(ellps='WGS84').inv(*LL_LonLat, *UL_LonLat)[0]
forward_az_right = pyproj.Geod(ellps='WGS84').inv(*LR_LonLat, *UR_LonLat)[0]
......@@ -1189,9 +1184,10 @@ def calc_VAA_using_fullSceneCornerLonLat(fullSceneCornerLonLat, orbit_params):
def calc_VZA_array(shape_fullArr, arr_pos, fullSceneCornerPos, viewing_angle, FOV, logger, meshwidth=1,
nodata_mask=None, outFill=None):
"""Calculates viewing zenith angles for each pixel in the dataset by solving
an equation system with 4 variables for each image corner:
VZA = a + b*col + c*row + d*col*row
"""Calculate viewing zenith angles for each pixel in the dataset.
By solving an equation system with 4 variables for each image corner: VZA = a + b*col + c*row + d*col*row.
:param shape_fullArr:
:param arr_pos:
:param fullSceneCornerPos:
......@@ -1203,7 +1199,6 @@ def calc_VZA_array(shape_fullArr, arr_pos, fullSceneCornerPos, viewing_angle, FO
:param nodata_mask: <numpy array>, used for declaring nodata values in the output VZA array
:param outFill: the value that is assigned to nodata area in the output VZA array
"""
# FIXME in case of Sentinel-2 the viewing_angle corresponds to the center point of the image footprint
# FIXME (trueDataCornerPos)
# FIXME => the algorithm must use the center viewing angle + orbit inclination and must calculate the FOV to be used
......@@ -1283,11 +1278,11 @@ def calc_AcqTime_array(shape_fullArr, arr_pos, AcqDate, CenterAcqTime, fullScene
def calc_SZA_SAA(date, lon, lat): # not used anymore since pyorbital is more precise
"""Calculates solar zenith and azimuth angle using pyephem.
:param date:
:param lon:
:param lat:
"""
obsv = ephem.Observer()
obsv.lon, obsv.lat = str(lon), str(lat)
obsv.date = date
......@@ -1303,6 +1298,7 @@ def calc_SZA_SAA_array(shape_fullArr, arr_pos, AcqDate, CenterAcqTime, fullScene
overpassDurationSec, logger, meshwidth=1, nodata_mask=None, outFill=None, accurracy='coarse',
lonlat_arr=None):
"""Calculates solar zenith and azimuth angles for each pixel in the dataset using pyorbital.
:param shape_fullArr:
:param arr_pos:
:param AcqDate:
......@@ -1323,7 +1319,6 @@ def calc_SZA_SAA_array(shape_fullArr, arr_pos, AcqDate, CenterAcqTime, fullScene
SZA/SAA = a + b*col + c*row + d*col*row.
:param lonlat_arr:
"""
if nodata_mask is not None:
assert isinstance(nodata_mask, (GeoArray, np.ndarray)), \
"'nodata_mask' must be a numpy array or an instance of GeoArray. Got %s" % type(nodata_mask)
......@@ -1402,7 +1397,6 @@ def calc_RAA_array(SAA_array, VAA_array, nodata_mask=None, outFill=None):
:param outFill: the value to be used to fill areas outside the actual image bounds
:return:
"""
if nodata_mask is not None:
assert isinstance(nodata_mask, (GeoArray, np.ndarray)), \
"'nodata_mask' must be a numpy array or an instance of GeoArray. Got %s" % type(nodata_mask)
......@@ -1476,7 +1470,6 @@ def clip_array_using_mapBounds(array, bounds, im_prj, im_gt, fillVal=0):
:param im_gt:
:param fillVal:
"""
print(bounds)
# get target bounds on the same grid like the input array
tgt_xmin, tgt_ymin, tgt_xmax, tgt_ymax = snap_bounds_to_pixGrid(bounds, im_gt)
......
......@@ -351,16 +351,12 @@ class Usecase:
self.target_FWHM = query_vir('band_width', self.virtual_sensor_id)
self.target_gsd = query_vir('spatial_resolution',
self.virtual_sensor_id) # table features only 1 value for X/Y-dims
self.target_gsd = \
[self.target_gsd, self.target_gsd] if type(self.target_gsd) in [int, float] else self.target_gsd
self.target_gsd = xgsd, ygsd = \
[self.target_gsd]*2 if isinstance(self.target_gsd, (int, float)) else self.target_gsd
self.EPSG = query_vir('projection_epsg', self.virtual_sensor_id)
self.spatial_ref_gridx = np.array(np.arange(self.target_gsd[0] / 2.,
self.target_gsd[0] / 2. + 2 * self.target_gsd[0],
self.target_gsd[0])) # e.g. [15,45
self.spatial_ref_gridy = np.array(np.arange(self.target_gsd[1] / 2.,
self.target_gsd[1] / 2. + 2 * self.target_gsd[1],
self.target_gsd[1]))
self.spatial_ref_gridx = np.arange(xgsd / 2., xgsd / 2. + 2 * xgsd, xgsd) # e.g. [15, 45]
self.spatial_ref_gridy = np.arange(ygsd / 2., ygsd / 2. + 2 * ygsd, ygsd)
self.scale_factor_TOARef = int(query_cfg('scale_factor_TOARef'))
self.scale_factor_BOARef = int(query_cfg('scale_factor_BOARef'))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment