Commit 9e2baf50 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Fixed bug regarding matplotlib backend. PEP-8 editing.

Former-commit-id: 5107ea13
Former-commit-id: b3925713
parent e2d189ef
# -*- 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,10 +11,14 @@ from . import processing
from . import config
from .processing.process_controller import process_controller
__all__ = ['algorithms',
'io',
'misc',
'processing',
'config',
'process_controller',
]
__author__ = """Daniel Scheffler"""
__email__ = 'daniel.scheffler@gfz-potsdam.de'
__version__ = '0.6.3'
__versionalias__ = '20170915.01'
__all__ = ['algorithms',
'io',
'misc',
'processing',
'config',
'process_controller',
]
......@@ -46,27 +46,27 @@ except ImportError:
from gdalconst import *
import pyproj
from pyorbital import astronomy
from pyorbital import astronomy
import ephem
from shapely.geometry import MultiPoint, Polygon
from shapely.ops import cascaded_union
from geoarray import GeoArray
from py_tools_ds.geo.coord_grid import snap_bounds_to_pixGrid
from py_tools_ds.geo.coord_trafo import transform_utm_to_wgs84, transform_wgs84_to_utm, mapXY2imXY, imXY2mapXY
from py_tools_ds.geo.projection import get_UTMzone, EPSG2WKT, isProjectedOrGeographic
from py_tools_ds.geo.coord_grid import snap_bounds_to_pixGrid
from py_tools_ds.geo.coord_trafo import transform_utm_to_wgs84, transform_wgs84_to_utm, mapXY2imXY, imXY2mapXY
from py_tools_ds.geo.projection import get_UTMzone, EPSG2WKT, isProjectedOrGeographic
from py_tools_ds.geo.raster.reproject import warp_ndarray
from ..config import GMS_config as CFG
from ..misc import helper_functions as HLP_F
from ..io import envifilehandling as ef
from ..misc import helper_functions as HLP_F
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,20 +104,21 @@ 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()
geodata = None
self.filepath, self.filename = None, self.inDs.GetDescription()
self.shortname, self.extension = os.path.splitext(self.filename)
# ****OBJECT ATTRIBUTES***************************************************
self.workspace = "/dev/shm/GeoMultiSens/GEOPROCESSING_temp/" if workspace is None else workspace
self.workspace = "/dev/shm/GeoMultiSens/GEOPROCESSING_temp/" if workspace is None else workspace
if v:
self.logger.debug("\n--")
self.logger.debug("\ttemporary geoprocessing workspace", self.workspace)
......@@ -125,16 +126,16 @@ class GEOPROCESSING(object):
# --Getting columns, rows and number of bands of inputdataset
self.desc = self.inDs.GetDescription()
if self.subset is None or self.subset[0] == 'cube':
self.cols = self.inDs.RasterXSize
self.rows = self.inDs.RasterYSize
self.bands = self.inDs.RasterCount
self.colStart = 0
self.rowStart = 0
self.cols = self.inDs.RasterXSize
self.rows = self.inDs.RasterYSize
self.bands = self.inDs.RasterCount
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
......@@ -146,28 +147,28 @@ class GEOPROCESSING(object):
%s\n\tbands[bands]: %s\n\tDescription[desc]: %s""" % (self.cols, self.rows, self.bands, self.desc))
# --Getting driver infos of inputdataset
self.drname_s = self.inDs.GetDriver().ShortName
self.drname_l = self.inDs.GetDriver().LongName
self.drhelp = self.inDs.GetDriver().HelpTopic
self.drname_s = self.inDs.GetDriver().ShortName
self.drname_l = self.inDs.GetDriver().LongName
self.drhelp = self.inDs.GetDriver().HelpTopic
self.DriverProp = [self.drname_s, self.drname_l, self.drhelp]
if v is not None:
self.logger.info(
"\tDriverShortName[drname_s]: %s\n\tDriverLongName[drnam_l]: %s\n\tDriverHelpWebsite[drhelp]: %s" % (
self.drname_s, self.drname_l, self.drhelp))
self.drname_s, self.drname_l, self.drhelp))
# --Getting Georeference info of inputdataset
self.geotransform = self.inDs.GetGeoTransform()
self.projection = self.inDs.GetProjection()
self.originX = self.geotransform[0]
self.originY = self.geotransform[3]
self.pixelWidth = self.geotransform[1]
self.pixelHeight = self.geotransform[5]
self.rot1 = self.geotransform[2]
self.rot2 = self.geotransform[4]
self.extent = [self.originX, self.originY, self.originX + (self.cols * self.pixelWidth),
self.originY + (self.rows * self.pixelHeight)] # [ulx, uly, lrx, lry]
self.GeoProp = [self.originX, self.originY, self.pixelWidth, self.pixelHeight, self.rot1, self.rot2]
self.projection = self.inDs.GetProjection()
self.originX = self.geotransform[0]
self.originY = self.geotransform[3]
self.pixelWidth = self.geotransform[1]
self.pixelHeight = self.geotransform[5]
self.rot1 = self.geotransform[2]
self.rot2 = self.geotransform[4]
self.extent = [self.originX, self.originY, self.originX + (self.cols * self.pixelWidth),
self.originY + (self.rows * self.pixelHeight)] # [ulx, uly, lrx, lry]
self.GeoProp = [self.originX, self.originY, self.pixelWidth, self.pixelHeight, self.rot1, self.rot2]
if v is not None:
self.logger.info("\toriginX[originX]:", self.originX, "\n\toriginY[originY]:", self.originY,
......@@ -178,48 +179,49 @@ 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()
self.filepath, self.filename = os.path.split(self.desc)
self.desc = self.inDs.GetDescription()
self.filepath, self.filename = os.path.split(self.desc)
self.shortname, self.extension = os.path.splitext(self.filename)
if self.subset is None or self.subset == 'cube':
self.cols = self.inDs.RasterXSize
self.rows = self.inDs.RasterYSize
self.bands = self.inDs.RasterCount
self.colStart = 0
self.rowStart = 0
self.cols = self.inDs.RasterXSize
self.rows = self.inDs.RasterYSize
self.bands = self.inDs.RasterCount
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()
self.DataProp = [self.cols, self.rows, self.bands, self.colStart, self.rowStart, self.bandStart]
self.drname_s = self.inDs.GetDriver().ShortName
self.drname_l = self.inDs.GetDriver().LongName
self.drhelp = self.inDs.GetDriver().HelpTopic
self.DataProp = [self.cols, self.rows, self.bands, self.colStart, self.rowStart, self.bandStart]
self.drname_s = self.inDs.GetDriver().ShortName
self.drname_l = self.inDs.GetDriver().LongName
self.drhelp = self.inDs.GetDriver().HelpTopic
self.DriverProp = [self.drname_s, self.drname_l, self.drhelp]
self.geotransform = self.inDs.GetGeoTransform()
self.projection = self.inDs.GetProjection()
self.originX = self.geotransform[0]
self.originY = self.geotransform[3]
self.pixelWidth = self.geotransform[1]
self.pixelHeight = self.geotransform[5]
self.rot1 = self.geotransform[2] # FIXME check
self.rot2 = self.geotransform[4]
self.extent = [self.originX, self.originY, self.originX + (self.cols * self.pixelWidth),
self.originY + (self.rows * self.pixelHeight)] # [ulx, uly, lrx, lry]
self.GeoProp = [self.originX, self.originY, self.pixelWidth, self.pixelHeight, self.rot1, self.rot2]
self.projection = self.inDs.GetProjection()
self.originX = self.geotransform[0]
self.originY = self.geotransform[3]
self.pixelWidth = self.geotransform[1]
self.pixelHeight = self.geotransform[5]
self.rot1 = self.geotransform[2] # FIXME check
self.rot2 = self.geotransform[4]
self.extent = [self.originX, self.originY, self.originX + (self.cols * self.pixelWidth),
self.originY + (self.rows * self.pixelHeight)] # [ulx, uly, lrx, lry]
self.GeoProp = [self.originX, self.originY, self.pixelWidth, self.pixelHeight, self.rot1, self.rot2]
self.Prop = {'DataProp': self.DataProp, 'DriverProp': self.DriverProp, 'GeoProp': self.GeoProp}
......@@ -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!
......@@ -491,7 +495,7 @@ class GEOPROCESSING(object):
self.logger.info("\t##\n\tband" + str(x + 1) + " of " + str(self.bands) + " completed")
# --Close Gdalobjects python: GdalObject = None -> behaves like GDALClose() in C
band = None
band = None
outBand = None
# self.inDs = None
......@@ -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
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]
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]
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,
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)
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)
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)
......@@ -657,7 +663,7 @@ class GEOPROCESSING(object):
mask = np.greater(data_Red + data_NIR, 0)
# exclude not interpretable values (for example saturated pixel)
if exclValues is None:
mask = mask # mask = mask_
mask = mask # mask = mask_
else:
# mask = np.where(np.any([data_Red == np.all(exclValues),data_NIR == np.all(exclValues)]),0,mask_)
x = 0
......@@ -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,10 +930,13 @@ 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]
# distance in [m]/cols
self.logger.info("While adding projection the ground sampling distance had to be calculated. It has "
"been set to %s meters." % gResolution)
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 "
"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]
self.inDs.SetGeoTransform(GT)
......@@ -968,60 +976,59 @@ 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)
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