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

Fixed bug regarding matplotlib backend. PEP-8 editing.


Former-commit-id: 5107ea13
parent 8c82beff
# -*- coding: utf-8 -*-
__author__ = """Daniel Scheffler"""
__email__ = 'daniel.scheffler@gfz-potsdam.de'
__version__ = '0.6.2'
__versionalias__ = '20170906.01'
import os
if 'MPLBACKEND' not in os.environ:
os.environ['MPLBACKEND'] = 'Agg'
from . import algorithms
from . import io
......@@ -13,6 +11,10 @@ from . import processing
from . import config
from .processing.process_controller import process_controller
__author__ = """Daniel Scheffler"""
__email__ = 'daniel.scheffler@gfz-potsdam.de'
__version__ = '0.6.3'
__versionalias__ = '20170915.01'
__all__ = ['algorithms',
'io',
'misc',
......
......@@ -63,10 +63,10 @@ from ..io import envifilehandling as ef
from ..misc.definition_dicts import get_outFillZeroSaturated
class GEOPROCESSING(object):
"""****CREATE OBJECT ****************************************************"""
"""****OBJECT CONSTRUCTOR**************************************************"""
def __init__(self, geodata, logger, workspace=None, subset=None, v=None):
self.logger = logger
self.subset = subset
......@@ -90,9 +90,9 @@ class GEOPROCESSING(object):
gdal.AllRegister()
assert type(geodata) in [str,gdal.Dataset], \
"GEOP: The argument 'geodata' has to be a path or a gdal.Dataset. Got %s." %type(geodata)
if isinstance(geodata,str):
assert type(geodata) in [str, gdal.Dataset], \
"GEOP: The argument 'geodata' has to be a path or a gdal.Dataset. Got %s." % type(geodata)
if isinstance(geodata, str):
# ----open files Open(<filename>, <GDALAccess>)
# open raster to gdal object
self.inDs = (gdal.Open(geodata, GA_ReadOnly))
......@@ -104,13 +104,14 @@ class GEOPROCESSING(object):
if not geodata.startswith('/vsi'):
path2check = geodata
else:
gdal_prefix_dict = {'/vsizip':'.zip', '/vsitar':'.tar', '/vsitar':'.tar.gz', '/vsitar' '.gz':'/vsigzip'}
p1 = [geodata.split(i)[0]+i for i in ['.zip','.tar','.tar.gz','.gz','.tgz'] \
gdal_prefix_dict = {'/vsizip': '.zip', '/vsitar': '.tar', '/vsitar': '.tar.gz',
'/vsitar' '.gz': '/vsigzip'}
p1 = [geodata.split(i)[0] + i for i in ['.zip', '.tar', '.tar.gz', '.gz', '.tgz'] \
if len(geodata.split(i)) > 1 and geodata.split(i)[1].startswith('/')][0]
path2check = os.path.abspath('.'+re.search('/vsi[\s\S]*(/[\s\S,.]*)',p1,re.I).group(1))
assert os.path.exists(path2check), "ERROR: data %s does not exist!" %path2check
assert self.inDs is not None, "ERROR: Could not open %s!" %self.filename
elif isinstance(geodata,gdal.Dataset):
path2check = os.path.abspath('.' + re.search('/vsi[\s\S]*(/[\s\S,.]*)', p1, re.I).group(1))
assert os.path.exists(path2check), "ERROR: data %s does not exist!" % path2check
assert self.inDs is not None, "ERROR: Could not open %s!" % self.filename
elif isinstance(geodata, gdal.Dataset):
self.inDs = geodata
geodata = None
self.filepath, self.filename = None, self.inDs.GetDescription()
......@@ -131,10 +132,10 @@ class GEOPROCESSING(object):
self.colStart = 0
self.rowStart = 0
self.bandStart = 0
self.colEnd = self.cols-1
self.rowEnd = self.rows-1
self.bandEnd = self.bands-1
self.bandsList = range(self.bandStart,self.bandEnd+1)
self.colEnd = self.cols - 1
self.rowEnd = self.rows - 1
self.bandEnd = self.bands - 1
self.bandsList = range(self.bandStart, self.bandEnd + 1)
else:
assert isinstance(self.subset, list) and len(self.subset) == 2, \
"Unable to initialize GEOP object due to bad arguments at subset keyword. Got %s." % self.subset
......@@ -178,11 +179,12 @@ class GEOPROCESSING(object):
self.Prop = {'DataProp': self.DataProp, 'DriverProp': self.DriverProp, 'GeoProp': self.GeoProp}
"""****OBJECT METHODS******************************************************"""
def subset_kwargs_to_cols_rows_bands_colStart_rowStart_bandStart(self):
shape_fullArr = [self.inDs.RasterYSize,self.inDs.RasterXSize,self.inDs.RasterCount]
self.rows,self.cols,self.bands,self.rowStart,self.rowEnd, self.colStart,self.colEnd,self.bandStart,\
self.bandEnd, self.bandsList = get_subsetProps_from_subsetArg(shape_fullArr,self.subset).values()
self.subset = self.subset if self.subset[0]!='cube' else None
shape_fullArr = [self.inDs.RasterYSize, self.inDs.RasterXSize, self.inDs.RasterCount]
self.rows, self.cols, self.bands, self.rowStart, self.rowEnd, self.colStart, self.colEnd, self.bandStart, \
self.bandEnd, self.bandsList = get_subsetProps_from_subsetArg(shape_fullArr, self.subset).values()
self.subset = self.subset if self.subset[0] != 'cube' else None
def update_dataset_related_attributes(self):
self.desc = self.inDs.GetDescription()
......@@ -196,10 +198,10 @@ class GEOPROCESSING(object):
self.colStart = 0
self.rowStart = 0
self.bandStart = 0
self.colEnd = self.cols-1
self.rowEnd = self.rows-1
self.bandEnd = self.bands-1
self.bandsList = range(self.bandStart,self.bandEnd+1)
self.colEnd = self.cols - 1
self.rowEnd = self.rows - 1
self.bandEnd = self.bands - 1
self.bandsList = range(self.bandStart, self.bandEnd + 1)
else:
self.subset_kwargs_to_cols_rows_bands_colStart_rowStart_bandStart()
......@@ -247,7 +249,7 @@ class GEOPROCESSING(object):
"""
# Wenn Zusatzfunktionen fertig dann hier einfach den Dataset einmal durchlaufen lassen und Properties in
# Liste abspeichern. Vllt auch noch PrintHistogramm methode???"
self.BandProps = [GetBandProps(self.inDs.GetRasterBand(x+1)) for x in self.bandsList]
self.BandProps = [GetBandProps(self.inDs.GetRasterBand(x + 1)) for x in self.bandsList]
def DN2Rad(self, offsets, gains, inFill=None, inZero=None, inSaturated=None, cutNeg=True):
""" ----METHOD_3----------------------------------------------------------
......@@ -267,20 +269,22 @@ class GEOPROCESSING(object):
self.conversion_type_optical = 'Rad'
tempArr = gdalnumeric.LoadFile(self.desc) if not self.subset or self.subset[0] == 'cube' else \
gdalnumeric.LoadFile(self.desc,self.colStart,self.rowStart,self.cols,self.rows)
gdalnumeric.LoadFile(self.desc, self.colStart, self.rowStart, self.cols, self.rows)
assert tempArr is not None, 'Error reading file %s.' %self.desc
assert tempArr is not None, 'Error reading file %s.' % self.desc
# select needed offsets according to self.bandsList (needed if optical and thermal data are processed at once)
offsets, gains = [list(np.array(par)[self.bandsList]) for par in [offsets, gains]]
inArray = np.swapaxes(np.swapaxes(tempArr,0,2),0,1)
inArray = np.take(inArray,self.bandsList,axis=2) if self.subset and self.subset[0]=='custom' else inArray # FIXME ineffektiv, da unnötige Bänder eingelesen werden
rad = DN2Rad(inArray,offsets, gains, inFill=inFill, inZero=inZero, inSaturated=inSaturated,cutNeg=cutNeg)
rad = np.swapaxes(np.swapaxes(rad,0,2),1,2) # [rows,cols,bands] => [bands,rows,cols]
inArray = np.swapaxes(np.swapaxes(tempArr, 0, 2), 0, 1)
inArray = np.take(inArray, self.bandsList, axis=2) if self.subset and self.subset[
0] == 'custom' else inArray # FIXME ineffektiv, da unnötige Bänder eingelesen werden
rad = DN2Rad(inArray, offsets, gains, inFill=inFill, inZero=inZero, inSaturated=inSaturated, cutNeg=cutNeg)
rad = np.swapaxes(np.swapaxes(rad, 0, 2), 1, 2) # [rows,cols,bands] => [bands,rows,cols]
filename_in_MEM = os.path.splitext(self.filename)[0] + '__' + (str(self.subset[1]) if self.subset
and self.subset[0] != 'cube' else '') + 'rad.mem'
self.inDs = gdalnumeric.SaveArray(rad,filename_in_MEM,format='MEM',prototype=self.desc)
and self.subset[
0] != 'cube' else '') + 'rad.mem'
self.inDs = gdalnumeric.SaveArray(rad, filename_in_MEM, format='MEM', prototype=self.desc)
self.update_dataset_related_attributes()
def DN2TOARef(self, offsets, gains, irradiance, zenith, earthSunDist, inFill=None, inZero=None, inSaturated=None,
......@@ -306,23 +310,24 @@ class GEOPROCESSING(object):
self.conversion_type_optical = 'TOA_Ref'
tempArr = gdalnumeric.LoadFile(self.desc) if not self.subset or self.subset[0] == 'cube' else \
gdalnumeric.LoadFile(self.desc,self.colStart,self.rowStart,self.cols,self.rows)
assert tempArr is not None, 'Error reading file %s.' %self.desc
gdalnumeric.LoadFile(self.desc, self.colStart, self.rowStart, self.cols, self.rows)
assert tempArr is not None, 'Error reading file %s.' % self.desc
# select needed offsets according to self.bandsList (needed if optical and thermal data are processed at once)
offsets, gains, irradiance = [list(np.array(par)[self.bandsList]) for par in [offsets, gains, irradiance]]
inArray = np.swapaxes(np.swapaxes(tempArr,0,2),0,1)
inArray = np.take(inArray,self.bandsList,axis=2) if self.subset and self.subset[0]=='custom' else inArray # FIXME ineffektiv, da unnötige Bänder eingelesen werden
TOA_ref = DN2TOARef(inArray,offsets, gains, irradiance, zenith, earthSunDist,
inFill=inFill, inZero=inZero, inSaturated=inSaturated,cutNeg=cutNeg)
TOA_ref = np.swapaxes(np.swapaxes(TOA_ref,0,2),1,2) # [rows,cols,bands] => [bands,rows,cols]
inArray = np.swapaxes(np.swapaxes(tempArr, 0, 2), 0, 1)
inArray = np.take(inArray, self.bandsList, axis=2) if self.subset and self.subset[
0] == 'custom' else inArray # FIXME ineffektiv, da unnötige Bänder eingelesen werden
TOA_ref = DN2TOARef(inArray, offsets, gains, irradiance, zenith, earthSunDist,
inFill=inFill, inZero=inZero, inSaturated=inSaturated, cutNeg=cutNeg)
TOA_ref = np.swapaxes(np.swapaxes(TOA_ref, 0, 2), 1, 2) # [rows,cols,bands] => [bands,rows,cols]
filename_in_MEM = os.path.splitext(self.filename)[0] + '__' + (str(self.subset[1]) if self.subset
and self.subset[0] != 'cube' else '') + '_TOA_ref.mem'
self.inDs = gdalnumeric.SaveArray(TOA_ref,filename_in_MEM,format='MEM',prototype=self.desc)
and self.subset[
0] != 'cube' else '') + '_TOA_ref.mem'
self.inDs = gdalnumeric.SaveArray(TOA_ref, filename_in_MEM, format='MEM', prototype=self.desc)
self.update_dataset_related_attributes()
def DN2TOARefOLI(self, offsetsRef, gainsRef, zenith, outPath=None, fill=None, zero=None, saturated=None,
cutNeg=True, optshift=None, v=0):
"""----METHOD_3a----------------------------------------------------------
......@@ -423,7 +428,6 @@ class GEOPROCESSING(object):
# shift of the UL-coordinates about a half pixel to the northwest
outDs.SetProjection(self.inDs.GetProjection())
# --Convert to TOA Reflectance
for x, out_idx in zip(self.bandsList, range(self.bands)):
# self.bandsList == range(self.bands) if subset != custom!
......@@ -510,20 +514,22 @@ class GEOPROCESSING(object):
self.conversion_type_optical = 'Temp'
tempArr = gdalnumeric.LoadFile(self.desc) if not self.subset or self.subset[0] == 'cube' else \
gdalnumeric.LoadFile(self.desc,self.colStart,self.rowStart,self.cols,self.rows)
assert tempArr is not None, 'Error reading file %s.' %self.desc
gdalnumeric.LoadFile(self.desc, self.colStart, self.rowStart, self.cols, self.rows)
assert tempArr is not None, 'Error reading file %s.' % self.desc
# select needed offsets according to self.bandsList (needed if optical and thermal data are processed at once)
K1, K2 = [list(np.array(par)[self.bandsList]) for par in [K1, K2]]
inArray = np.swapaxes(np.swapaxes(tempArr,0,2),0,1)
inArray = np.take(inArray,self.bandsList,axis=2) if self.subset and self.subset[0]=='custom' else inArray # FIXME ineffektiv, da unnötige Bänder eingelesen werden
inArray = np.swapaxes(np.swapaxes(tempArr, 0, 2), 0, 1)
inArray = np.take(inArray, self.bandsList, axis=2) if self.subset and self.subset[
0] == 'custom' else inArray # FIXME ineffektiv, da unnötige Bänder eingelesen werden
Kelvin = TOARad2Kelvin_fastforward(inArray, K1, K2, emissivity=emissivity,
inFill=inFill, inZero=inZero, inSaturated=inSaturated)
Kelvin = np.swapaxes(np.swapaxes(Kelvin,0,2),1,2) # [rows,cols,bands] => [bands,rows,cols]
Kelvin = np.swapaxes(np.swapaxes(Kelvin, 0, 2), 1, 2) # [rows,cols,bands] => [bands,rows,cols]
filename_in_MEM = os.path.splitext(self.filename)[0] + '__' + (str(self.subset[1]) if self.subset
and self.subset[0] != 'cube' else '') + '_Kelvin.mem'
self.inDs = gdalnumeric.SaveArray(Kelvin,filename_in_MEM,format='MEM',prototype=self.desc)
and self.subset[
0] != 'cube' else '') + '_Kelvin.mem'
self.inDs = gdalnumeric.SaveArray(Kelvin, filename_in_MEM, format='MEM', prototype=self.desc)
self.update_dataset_related_attributes()
def DN2DegreesCelsius_fastforward(self, offsets, gains, K1, K2, emissivity=0.95, inFill=None, inZero=None,
......@@ -543,20 +549,22 @@ class GEOPROCESSING(object):
self.conversion_type_optical = 'Temp'
tempArr = gdalnumeric.LoadFile(self.desc) if not self.subset or self.subset[0] == 'cube' else \
gdalnumeric.LoadFile(self.desc,self.colStart,self.rowStart,self.cols,self.rows)
assert tempArr is not None, 'Error reading file %s.' %self.desc
gdalnumeric.LoadFile(self.desc, self.colStart, self.rowStart, self.cols, self.rows)
assert tempArr is not None, 'Error reading file %s.' % self.desc
# select needed offsets according to self.bandsList (needed if optical and thermal data are processed at once)
offsets, gains, K1, K2 = [list(np.array(par)[self.bandsList]) for par in [offsets, gains, K1, K2]]
inArray = np.swapaxes(np.swapaxes(tempArr,0,2),0,1)
inArray = np.take(inArray,self.bandsList,axis=2) if self.subset and self.subset[0]=='custom' else inArray # FIXME ineffektiv, da unnötige Bänder eingelesen werden
degreesCelsius = DN2DegreesCelsius_fastforward(inArray,offsets, gains, K1, K2, emissivity=emissivity,
inArray = np.swapaxes(np.swapaxes(tempArr, 0, 2), 0, 1)
inArray = np.take(inArray, self.bandsList, axis=2) if self.subset and self.subset[
0] == 'custom' else inArray # FIXME ineffektiv, da unnötige Bänder eingelesen werden
degreesCelsius = DN2DegreesCelsius_fastforward(inArray, offsets, gains, K1, K2, emissivity=emissivity,
inFill=inFill, inZero=inZero, inSaturated=inSaturated)
degreesCelsius = np.swapaxes(np.swapaxes(degreesCelsius,0,2),1,2) # [rows,cols,bands] => [bands,rows,cols]
filename_in_MEM= os.path.splitext(self.filename)[0] + '__' + (str(self.subset[1]) if self.subset
and self.subset[0] != 'cube' else '') + '_degreesCelsius.mem'
self.inDs = gdalnumeric.SaveArray(degreesCelsius,filename_in_MEM,format='MEM',prototype=self.desc)
degreesCelsius = np.swapaxes(np.swapaxes(degreesCelsius, 0, 2), 1, 2) # [rows,cols,bands] => [bands,rows,cols]
filename_in_MEM = os.path.splitext(self.filename)[0] + '__' + (str(self.subset[1]) if self.subset
and self.subset[
0] != 'cube' else '') + '_degreesCelsius.mem'
self.inDs = gdalnumeric.SaveArray(degreesCelsius, filename_in_MEM, format='MEM', prototype=self.desc)
self.update_dataset_related_attributes()
def NDVI(self, band1, band2, outGDT="GDT_Byte", outPath=None, exclValues=None, v=0):
......@@ -614,7 +622,6 @@ class GEOPROCESSING(object):
print("\t %s" % outPath)
print("\t*exclValue: %s" % exclValues)
# --Create Outpufile
# Define name of outputfile
if outPath == self.workspace:
......@@ -641,7 +648,6 @@ class GEOPROCESSING(object):
outDs.SetGeoTransform(self.inDs.GetGeoTransform())
outDs.SetProjection(self.inDs.GetProjection())
# --Calculate NDVI
if self.subset is None or self.subset[0] == 'cube':
data_NIR = self.inDs.GetRasterBand(band1).ReadAsArray(0, 0, self.cols, self.rows).astype(np.float32)
......@@ -692,7 +698,6 @@ class GEOPROCESSING(object):
print("\t##\n\tNDVI calculation completed")
# --Close Gdalobjects python: GdalObject = None -> behaves like GDALClose() in C
band = None
outBand = None
......@@ -751,7 +756,8 @@ class GEOPROCESSING(object):
else: # max extent without background of ascending images
for i, row in enumerate(bandSum[::bin]):
line = np.nonzero(row > 0) # gibt die indices an, an denen überall Werte (Wert = >0) auf dieser Zeile stehen
line = np.nonzero(
row > 0) # gibt die indices an, an denen überall Werte (Wert = >0) auf dieser Zeile stehen
if len(line[0]) > 0:
# print "first image pixel of line %s: %s" %(i,line[0][0])
# print "last image pixel of line %s: %s" %(i,line[0][-1])
......@@ -819,8 +825,6 @@ class GEOPROCESSING(object):
if v == 1:
print("ulx,uly,lrx,lry:", ulx, uly, lrx, lry)
# --Create Outputfile
# Define name of outputfile
if outPath == self.workspace:
......@@ -902,7 +906,8 @@ class GEOPROCESSING(object):
elif CS_TYPE is not None and CS_DATUM is not None:
if CS_TYPE == 'UTM' and CS_DATUM == 'WGS84':
UTM_zone = int(1 + (CornerTieP_LonLat[0][0] + 180.0) / 6.0) if CS_UTM_ZONE is None else CS_UTM_ZONE
EPSG_code = int('326' + str(UTM_zone)) if CornerTieP_LonLat[0][1] > 0.0 else int('327' + str(UTM_zone))
EPSG_code = int('326' + str(UTM_zone)) if CornerTieP_LonLat[0][1] > 0.0 else int(
'327' + str(UTM_zone))
elif CS_TYPE == 'LonLat' and CS_DATUM == 'WGS84':
EPSG_code = 4326
else:
......@@ -925,9 +930,12 @@ class GEOPROCESSING(object):
"'gResolution' or 'shape_fullArr'!")
else:
gResolution = float(pyproj.Geod(ellps=CS_DATUM).inv(CornerTieP_LonLat[0][0],
CornerTieP_LonLat[0][1],CornerTieP_LonLat[1][0],CornerTieP_LonLat[1][1])[2]) / shape_fullArr[1]
CornerTieP_LonLat[0][1],
CornerTieP_LonLat[1][0],
CornerTieP_LonLat[1][1])[2]) / shape_fullArr[1]
# distance in [m]/cols
self.logger.info("While adding projection the ground sampling distance had to be calculated. It has "
self.logger.info(
"While adding projection the ground sampling distance had to be calculated. It has "
"been set to %s meters." % gResolution)
self.pixelWidth, self.pixelHeight, self.rot1, self.rot2 = (gResolution, gResolution * -1., 0, 0)
GT = [self.originX, self.pixelWidth, self.rot1, self.originY, self.rot2, self.pixelHeight]
......@@ -968,55 +976,54 @@ class GEOPROCESSING(object):
if 'GDAL_DATA' not in os.environ:
# Prevents "Unable to open EPSG support file gcs.csv".
self.logger.warning("GDAL_DATA variable not set. Setting it to '%s'."
%sys.executable.split('anaconda')[0] + 'anaconda/share/gdal')
anaconda_path = re.search('([\w+\S\s]*anaconda[\w+\S\s]*?)/',sys.executable,re.I).group(1)
os.environ['GDAL_DATA'] = os.path.join(anaconda_path,'share/gdal')
if not re.search('/anaconda[\w+\S\s]*?/share/gdal', os.environ['GDAL_DATA'],re.I):
% sys.executable.split('anaconda')[0] + 'anaconda/share/gdal')
anaconda_path = re.search('([\w+\S\s]*anaconda[\w+\S\s]*?)/', sys.executable, re.I).group(1)
os.environ['GDAL_DATA'] = os.path.join(anaconda_path, 'share/gdal')
if not re.search('/anaconda[\w+\S\s]*?/share/gdal', os.environ['GDAL_DATA'], re.I):
self.logger.warning('GDAL_DATA variable seems to be incorrectly set since Anaconda´s Python is '
'executed and GDAL_DATA = %s' %os.environ['GDAL_DATA'])
'executed and GDAL_DATA = %s' % os.environ['GDAL_DATA'])
assert mode in ['rasterio','GDAL'], "The 'mode' argument must be set to 'rasterio' or 'GDAL'."
assert mode in ['rasterio', 'GDAL'], "The 'mode' argument must be set to 'rasterio' or 'GDAL'."
if mode == 'rasterio':
"""needs no temporary files but does not support multiprocessing"""
out_prj = EPSG2WKT(dst_EPSG_code)
proj_geo = isProjectedOrGeographic(self.projection)
assert proj_geo in ['projected','geographic']
TieP = TieP if proj_geo=='geographic' else \
[transform_wgs84_to_utm(LonLat[0],LonLat[1],get_UTMzone(prj=self.projection)) for LonLat in TieP]
xmin,xmax,ymin,ymax = HLP_F.corner_coord_to_minmax(TieP)
t0= time.time()
in_arr = np.swapaxes(np.swapaxes(gdalnumeric.LoadFile(self.desc),0,2),0,1)
print('reading time', time.time()-t0)
assert proj_geo in ['projected', 'geographic']
TieP = TieP if proj_geo == 'geographic' else \
[transform_wgs84_to_utm(LonLat[0], LonLat[1], get_UTMzone(prj=self.projection)) for LonLat in TieP]
xmin, xmax, ymin, ymax = HLP_F.corner_coord_to_minmax(TieP)
t0 = time.time()
in_arr = np.swapaxes(np.swapaxes(gdalnumeric.LoadFile(self.desc), 0, 2), 0, 1)
print('reading time', time.time() - t0)
if inFill is None: inFill = get_outFillZeroSaturated(in_arr.dtype)[0]
out_nodataVal = get_outFillZeroSaturated(in_arr.dtype)[0]
t0= time.time()
t0 = time.time()
warped, gt, prj = warp_ndarray(in_arr, self.geotransform, self.projection, out_prj,
out_bounds=([xmin, ymin, xmax, ymax]), rspAlg='cubic',
in_nodata=inFill,out_nodata=out_nodataVal) [0]
in_nodata=inFill, out_nodata=out_nodataVal)[0]
# multiprocessing: not implemented further because multiproceesing within Python Mappers is not possible
#args = [(in_arr[:,:,i],self.geotransform,self.projection,out_prj,([xmin, ymin, xmax, ymax]),
# args = [(in_arr[:,:,i],self.geotransform,self.projection,out_prj,([xmin, ymin, xmax, ymax]),
# 2,in_nodataVal) for i in range(in_arr.shape[2])]
#import multiprocessing
#with multiprocessing.Pool() as pool:
# import multiprocessing
# with multiprocessing.Pool() as pool:
# results = pool.map(warp_mp,args)
print('warping time', time.time()-t0)
#from spectral.io import envi
#envi.save_image('/dev/shm/test_warped.hdr',warped,dtype=str(np.dtype(warped.dtype)),
print('warping time', time.time() - t0)
# from spectral.io import envi
# envi.save_image('/dev/shm/test_warped.hdr',warped,dtype=str(np.dtype(warped.dtype)),
# force=True,ext='.bsq',interleave='bsq')
else: # mode == 'GDAL'
"""needs temporary files but does support multiprocessing and configuring cache size"""
t0 = time.time()
in_dtype = gdalnumeric.LoadFile(self.desc,0,0,1,1).dtype
in_dtype = gdalnumeric.LoadFile(self.desc, 0, 0, 1, 1).dtype
if inFill is None: inFill = get_outFillZeroSaturated(in_dtype)[0]
out_nodataVal = get_outFillZeroSaturated(in_dtype)[0]
gcps = ' '.join(['-gcp %s %s %s %s' %tuple(gcp) for gcp in [gcp_ul,gcp_ur,gcp_ll,gcp_lr]])
gcps = ' '.join(['-gcp %s %s %s %s' % tuple(gcp) for gcp in [gcp_ul, gcp_ur, gcp_ll, gcp_lr]])
if use_workspace:
inFile = self.desc
......@@ -1031,7 +1038,7 @@ class GEOPROCESSING(object):
# %(dst_EPSG_code, in_nodataVal,out_nodataVal, translatedFile, warpedFile))
os.system('gdalwarp -of ENVI --config GDAL_CACHEMAX 2048 -wm 2048 -t_srs EPSG:%s -tps -r \
cubic -srcnodata %s -dstnodata %s -multi -overwrite -wo NUM_THREADS=%s -q %s %s' \
%(dst_EPSG_code,inFill,out_nodataVal,CFG.job.CPUs,translatedFile,warpedFile))
% (dst_EPSG_code, inFill, out_nodataVal, CFG.job.CPUs, translatedFile, warpedFile))
# import shutil
# shutil.copy(translatedFile, '//misc/hy5/scheffler/Skripte_Models/python/gms_preprocessing/sandbox/out/') ## only for bugfixing
# shutil.copy(translatedFile+'.hdr','//misc/hy5/scheffler/Skripte_Models/python/gms_preprocessing/sandbox/out/') ## only for bugfixing
......@@ -1054,12 +1061,11 @@ class GEOPROCESSING(object):
cubic -srcnodata %s -dstnodata %s -overwrite -multi -wo NUM_THREADS=%s -q %s %s' \
% (dst_EPSG_code, inFill, out_nodataVal, CFG.job.CPUs, translatedFile, warpedFile))
# print('warped')
print('GDAL warping time',time.time()-t0)
print('GDAL warping time', time.time() - t0)
self.inDs = gdal.OpenShared(warpedFile, GA_ReadOnly) \
if inFile.endswith('.vrt') else gdal.Open(warpedFile, GA_ReadOnly)
#print('read successful')
# print('read successful')
self.update_dataset_related_attributes()
# def get_row_column_bounds(self,arr_shape = None, arr_pos = None):
......@@ -1092,7 +1098,8 @@ class GEOPROCESSING(object):
=> 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)
UR_LonLat = self.pixelToWorldCoord_using_geotransform_and_projection((0, self.cols),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) # self.cols statt self.cols-1, um Koordinate am rechten Pixelrand zu berechnen
LL_LonLat = self.pixelToWorldCoord_using_geotransform_and_projection((self.rows, 0), targetProj)
LR_LonLat = self.pixelToWorldCoord_using_geotransform_and_projection((self.rows, self.cols), targetProj)
UL_LonLat = tuple([round(i, 10) for i in UL_LonLat])
......@@ -1101,18 +1108,17 @@ class GEOPROCESSING(object):
LR_LonLat = tuple([round(i, 10) for i in LR_LonLat])
return [UL_LonLat, UR_LonLat, LL_LonLat, LR_LonLat]
def calc_mask_data_nodata(self, array=None, custom_nodataVal = -9999):
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)
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)
def calc_mask_data_nodataOLD(self, custom_nodataVal = -9999):
def calc_mask_data_nodataOLD(self, custom_nodataVal=-9999):
"""Berechnet den Bildbereich, der von allen Kanälen abgedeckt wird.
:param custom_nodataVal:
"""
......@@ -1123,7 +1129,7 @@ class GEOPROCESSING(object):
# self.bandsList == range(self.bands) if subset != custom!
for x, out_idx in zip(self.bandsList, range(self.bands)):
band = self.inDs.GetRasterBand(x + 1) # Read inputband as numpyarray
cS,rS = [0,0] if self.subset is None or self.subset[0] == 'cube' else [self.colStart, self.rowStart]
cS, rS = [0, 0] if self.subset is None or self.subset[0] == 'cube' else [self.colStart, self.rowStart]
data = band.ReadAsArray(cS, rS, self.cols, self.rows).astype(np.int16)
# if self.bands > 1:
# img[:,:,out_idx] = data
......@@ -1156,7 +1162,7 @@ class GEOPROCESSING(object):
first_dataRow = list(rows_containing_data).index(True)
last_dataRow = self.rows - (list(reversed(rows_containing_data)).index(True) + 1)
if first_dataCol == 0 and first_dataRow == 0 and last_dataCol == self.cols-1 and last_dataRow == self.rows-1 \
if first_dataCol == 0 and first_dataRow == 0 and last_dataCol == self.cols - 1 and last_dataRow == self.rows - 1 \
and mask_1bit[0, 0] == True and mask_1bit[0, self.cols - 1] == True \
and mask_1bit[self.rows - 1, 0] == True and mask_1bit[self.rows - 1, self.cols - 1] == True:
UL, UR, LL, LR = (0, 0), (0, self.cols - 1), (self.rows - 1, 0), (self.rows - 1, self.cols - 1)
......@@ -1174,6 +1180,7 @@ class GEOPROCESSING(object):
[StartStopRows_in_first_dataCol, StartStopRows_in_last_dataCol, StartStopCols_in_first_dataRow,
StartStopCols_in_last_dataRow]]:
''' In case of cut image corners (e.g. ALOS AVNIR-2 Level-1B2): Calculation of fullSceneCornerPos outside of the image.'''
# print ('Image corners detected to be cut.')
def find_line_intersection_point(line1, line2):
xdiff = (line1[0][0] - line1[1][0], line2[0][0] - line2[1][0])
......@@ -1216,9 +1223,11 @@ class GEOPROCESSING(object):
distUL = [np.sqrt((0 - c[0]) ** 2 + (0 - c[1]) ** 2) for c in corners] # distance of each corner to UL
UL = corners[distUL.index(min(distUL))]
distUR = [np.sqrt((0 - c[0]) ** 2 + (self.cols - c[1]) ** 2) for c in corners] # distance of each corner to UR
distUR = [np.sqrt((0 - c[0]) ** 2 + (self.cols - c[1]) ** 2) for c in
corners] # distance of each corner to UR
UR = corners[distUR.index(min(distUR))]
distLL = [np.sqrt((self.rows - c[0]) ** 2 + (0 - c[1]) ** 2) for c in corners] # distance of each corner to LL
distLL = [np.sqrt((self.rows - c[0]) ** 2 + (0 - c[1]) ** 2) for c in
corners] # distance of each corner to LL
LL = corners[distLL.index(min(distLL))]
distLR = [np.sqrt((self.rows - c[0]) ** 2 + (self.cols - c[1]) ** 2) for c in
corners] # distance of each corner to LR
......@@ -1291,7 +1300,8 @@ class GEOPROCESSING(object):
DataType = band.DataType
GDALDataType = gdal.GetDataTypeName(DataType)
# dictionary to translate GDAL data types in corresponding numpy data types
dTypeDic = {"Byte": np.uint8, "UInt16": np.uint16, "Int16": np.int16, "UInt32": np.uint32, "Int32": np.int32,
dTypeDic = {"Byte": np.uint8, "UInt16": np.uint16, "Int16": np.int16, "UInt32": np.uint32,
"Int32": np.int32,
"Float32": np.float32, "Float64": np.float64}
else:
GDALDataType = dtype
......@@ -1416,6 +1426,7 @@ class GEOPROCESSING(object):
# data shifting
"""++++Conversion methods++++++++++++++++++++++++++++++++++++++++++++++++"""
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])
......@@ -1534,8 +1545,9 @@ class GEOPROCESSING(object):
print("\t extent:", extent)
if self.subset is not None and self.subset != 'cube': # rasObj instanced as subset
if [self.inDs.RasterYSize,self.inDs.RasterXSize] != [self.rows,self.cols]: # inDs does not represent subset
startpix = [startpix[0]+self.colStart, startpix[1]+self.rowStart]
if [self.inDs.RasterYSize, self.inDs.RasterXSize] != [self.rows,
self.cols]: # inDs does not represent subset
startpix = [startpix[0] + self.colStart, startpix[1] + self.rowStart]
extent = [self.cols, self.rows]
if self.subset is not None and self.subset[0] == 'custom' and \
......@@ -1547,13 +1559,14 @@ class GEOPROCESSING(object):
bands2process = range(self.bands)
np_dtype = convertGdalNumpyDataType(gdal.GetDataTypeName(self.inDs.GetRasterBand(1).DataType))
im = np.empty((self.rows,self.cols,len(bands2process)),np_dtype)
im = np.empty((self.rows, self.cols, len(bands2process)), np_dtype)
for out_idx, x in enumerate(bands2process):
im[:,:,out_idx] = self.inDs.GetRasterBand(x + 1).ReadAsArray(startpix[0], startpix[1], extent[0], extent[1])
im[:, :, out_idx] = self.inDs.GetRasterBand(x + 1).ReadAsArray(startpix[0], startpix[1], extent[0],
extent[1])
if direction == 1: # bands, rows, cols; so wie Theres es gebastelt hat
im = np.rollaxis(im,2)
im = np.rollaxis(im, 2)
elif direction == 2: # rows, bands, cols
im = np.swapaxes(im,1,2)
im = np.swapaxes(im, 1, 2)
elif direction == 3: # rows, cols, bands; so wie in PIL bzw. IDL benötigt
pass # bands, rows,cols
return im
......@@ -1622,7 +1635,8 @@ class GEOPROCESSING(object):
rotation1 = geotransform[2]
rotation2 = geotransform[4]
print("\toX:", originX, "\n\toY:", originY, "\n\tresolution: X:" + str(pixelWidth), "Y:" + str(pixelHeight),
print