Commit 76fe39f0 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

added new property 'mask_baddata' to GeoArray; some bugfixes and further...

added new property 'mask_baddata' to GeoArray; some bugfixes and further developments; new dtypes package

- added package 'dtypes' with submodule 'conversion' for performing data tyoe conversions

geo.raster.reproject:
- warp_ndarray(): input data types that are incompatible to GDAL are now automatically transformed to a GDAL compatible data type

io.raster.gdal:
- get_GDAL_ds_inmem(): added automatic data type conversion if not compatible to GDAL

io.raster.GeoArray:
- GeoArray:
    - __init__():
        - bugfix for rejecting subclasses of GeoArray
        - added attribute '_mask_baddata'
    - added property 'mask_baddata' + setter
- added class 'BadDataMask' (subclass of GeoArray)

- updated __version__
parent 42ae52b1
......@@ -15,7 +15,7 @@ __all__=[#'compatibility',
'similarity',
'GeoArray']
__version__ = '20161117_01'
__version__ = '20161122_01'
__author__='Daniel Scheffler'
# Validate GDAL version
......
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
from . import conversion
__all__=['conversion']
\ No newline at end of file
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
import numpy as np
try:
from osgeo import gdal
except ImportError:
import gdal
# dictionary to translate Numpy data types (strings) into corresponding GDAL data types,
# e.g. dTypeDic_NumPy2GDAL(str(np.dtype(np.uint8)))
dTypeDic_NumPy2GDAL = {'bool' : gdal.GDT_Byte,
'bool_' : gdal.GDT_Int32,
'int' : gdal.GDT_Int32,
'int8' : gdal.GDT_Int16,
'uint8' : gdal.GDT_Byte,
'uint16' : gdal.GDT_UInt16,
'int16' : gdal.GDT_Int16,
'uint32' : gdal.GDT_UInt32,
'int32' : gdal.GDT_Int32,
'int64' : gdal.GDT_Float64,
'float' : gdal.GDT_Float32,
'float16' : gdal.GDT_Float32,
'float32' : gdal.GDT_Float32,
'float64' : gdal.GDT_Float64
}
# dictionary to translate GDAL data types (strings) into corresponding numpy data types
dTypeDic_GDAL2Numpy = {gdal.GDT_Byte : np.uint8,
gdal.GDT_UInt16 : np.uint16,
gdal.GDT_Int16 : np.int16,
gdal.GDT_UInt32 : np.uint32,
gdal.GDT_Int32 : np.int32,
gdal.GDT_Float32 : np.float32,
gdal.GDT_Float64 : np.float64,
}
# dictionary to translate Numpy data types into GDAL compatible Numpy data types
dTypeDic_NumPy2GDALcompatible = \
dict(zip(dTypeDic_NumPy2GDAL.keys(),
[dTypeDic_GDAL2Numpy[dTypeDic_NumPy2GDAL[str(np.dtype(NDT))]] for NDT in dTypeDic_NumPy2GDAL.keys()]))
......@@ -18,6 +18,7 @@ from rasterio.warp import reproject as rio_reproject
from rasterio.warp import calculate_default_transform as rio_calc_transform
from rasterio.warp import Resampling
from ...dtypes.conversion import dTypeDic_NumPy2GDAL, dTypeDic_GDAL2Numpy
from ..projection import WKT2EPSG, isProjectedOrGeographic, prj_equal
from ..coord_trafo import pixelToLatLon
from ...io.raster.gdal import get_GDAL_ds_inmem
......@@ -25,16 +26,6 @@ from ...processing.progress_mon import ProgressBar
from ...compatibility.gdal import get_gdal_func
# dictionary to translate GDAL data types (strings) in corresponding numpy data types
dTypeDic_NumPy2GDAL = {'uint8' : gdal.GDT_Byte,
'uint16' : gdal.GDT_UInt16,
'int16' : gdal.GDT_Int16,
'uint32' : gdal.GDT_UInt32,
'int32' : gdal.GDT_Int32,
'float32' : gdal.GDT_Float32,
'float64' : gdal.GDT_Float64
}
def warp_ndarray_OLD(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None, outUL=None, out_res=None,
out_extent=None, out_dtype=None, rsp_alg=0, in_nodata=None, out_nodata=None, outExtent_within=True):
......@@ -278,10 +269,9 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=(
rspAlg='cubic'
get_SRS = lambda prjArg: prjArg if isinstance(prjArg,str) and prjArg.startswith('EPSG:') else \
'EPSG:%s'%prjArg if isinstance(prjArg,int) else \
prjArg
get_GDT = lambda DT: dTypeDic_NumPy2GDAL[str(np.dtype(DT))]
get_SRS = lambda prjArg: prjArg if isinstance(prjArg,str) and prjArg.startswith('EPSG:') else \
'EPSG:%s'%prjArg if isinstance(prjArg,int) else prjArg
get_GDT = lambda DT: dTypeDic_NumPy2GDAL[str(np.dtype(DT))]
# not yet implemented
cutlineDSName = 'data/cutline.vrt' #'/vsimem/cutline.shp' TODO cutline from OGR datasource. => implement input shapefile or Geopandas dataframe
......@@ -422,6 +412,10 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=(
# cleanup
in_ds = res_ds = None
# dtype check -> possibly dtype had to be changed for GDAL compatibility
if ndarray.dtype != res_arr.dtype:
res_arr = res_arr.astype(ndarray.dtype)
# test output
if out_prj and prj_equal(out_prj,4626):
assert -180 < res_gt[0] < 180 and -90 < res_gt[3] < 90, 'Testing of gdal_warp output failed.'
......
......@@ -41,9 +41,9 @@ from ...compatibility.gdal import get_gdal_func
def _alias_property(key):
return property(
lambda self: getattr(self, key),
lambda self: getattr(self, key),
lambda self, val: setattr(self, key, val),
lambda self: delattr(self, key))
lambda self: delattr(self, key))
class GeoArray(object):
......@@ -70,9 +70,13 @@ class GeoArray(object):
# FIXME implement compatibility to GDAL VRTs
if type(path_or_array) not in [str, np.ndarray, type(GeoArray)]:
raise ValueError("%s parameter 'arg' takes only string "
"or np.ndarray types. Got %s." %(self.__class__.__name__,type(path_or_array)))
if type(path_or_array) not in [str, np.ndarray, GeoArray] and not \
not issubclass(getattr(path_or_array,'__class__'), GeoArray):
raise ValueError("%s parameter 'arg' takes only string, np.ndarray or GeoArray types. Got %s."
%(self.__class__.__name__,type(path_or_array)))
if path_or_array is None:
raise ValueError("The %s parameter 'path_or_array' must not be None!" %self.__class__.__name__)
if isinstance(path_or_array, str):
assert ' ' not in path_or_array, "The given path contains whitespaces. This is not supported by GDAL."
......@@ -81,12 +85,12 @@ class GeoArray(object):
raise FileNotFoundError(path_or_array) if PY3 else FileNotFoundError_comp(path_or_array)
if type(path_or_array)==type(self):
if isinstance(path_or_array, GeoArray) or issubclass(getattr(path_or_array,'__class__'), GeoArray):
self.__dict__= path_or_array.__dict__
self._initParams = dict([x for x in locals().items() if x[0] != "self"])
self.geotransform = geotransform if geotransform else self.geotransform
self.projection = projection if projection else self.projection
self.bandnames = bandnames if bandnames else self.bandnames
self.bandnames = bandnames if bandnames else list(self.bandnames.values())
self._nodata = nodata if nodata is not None else self._nodata
self.progress = progress if progress else self.progress
self.q = q if q is not None else self.q
......@@ -106,6 +110,7 @@ class GeoArray(object):
self._dtype = None
self._nodata = nodata
self._mask_nodata = None
self._mask_baddata = None
self._footprint_poly = None
self._gdalDataset_meta_already_set = False
self._metadata = None
......@@ -143,6 +148,7 @@ class GeoArray(object):
@bandnames.setter
def bandnames(self, list_bandnames):
# type: (list)
assert len(list_bandnames) == self.bands, \
'Number of given bandnames does not match number of bands in array.'
assert len(list(set([type(b) for b in list_bandnames]))) == 1 and type(list_bandnames[0] == 'str'), \
......@@ -339,9 +345,44 @@ class GeoArray(object):
self._mask_nodata = maskarray.astype(np.uint8)
@property
def mask_baddata(self):
"""Returns the bad data mask for the associated image array if it has been explicitly previously. It can be set
by passing a file path, a numpy array or an instance of GeoArray to the setter of this property.
"""
return self._mask_baddata
@mask_baddata.setter
def mask_baddata(self, mask):
"""Set bad data mask.
:param mask: Can be a file path, a numpy array or an instance o GeoArray.
"""
if mask is not None:
geoArr_mask = BadDataMask(mask)
geoArr_mask.gt = geoArr_mask.gt if geoArr_mask.gt not in [None, [0, 1, 0, 0, 0, -1]] else self.gt
geoArr_mask.prj = geoArr_mask.prj if geoArr_mask.prj is None else self.prj
imName = "the %s '%s'" %(self.__class__.__name__, self.basename)
assert geoArr_mask.bands == 1, \
'Expected one single band as bad data mask for %s. Got %s bands.' % (self.basename, geoArr_mask.bands)
assert geoArr_mask.shape[:2] == self.shape[:2], 'The provided bad data mask must have the same number of ' \
'rows and columns as the %s itself.' %imName
assert geoArr_mask.gt == self.gt, \
'The geotransform of the given bad data mask for %s must match the geotransform of the %s itself. ' \
'Got %s.' %(imName, self.__class__.__name__, geoArr_mask.gt)
assert prj_equal(geoArr_mask.prj, self.prj), \
'The projection of the given bad data mask for the %s must match the projection of the %s itself.' \
%(imName, self.__class__.__name__)
self._mask_baddata = geoArr_mask
@property
def footprint_poly(self):
# FIXME should return polygon in image coordiates if no projection is available
# FIXME should return polygon in image coordinates if no projection is available
"""Get the footprint polygon of the associated image array (returns an instance of shapely.geometry.Polygon."""
if self._footprint_poly is not None:
......@@ -1075,6 +1116,33 @@ class GeoArray(object):
class BadDataMask(GeoArray):
def __init__(self, path_or_array, geotransform=None, projection=None, bandnames=None, nodata=None, progress=True,
q=False):
super(BadDataMask, self).__init__(path_or_array, geotransform=geotransform, projection=projection,
bandnames=bandnames, nodata=nodata, progress=progress, q=q)
self.arr = self.arr.astype(np.bool)
@property
def arr(self):
return self._arr
@arr.setter
def arr(self, ndarray):
assert isinstance(ndarray, np.ndarray), "'arr' can only be set to a numpy array!"
pixelVals_in_mask = sorted(list(np.unique(ndarray)))
assert len(pixelVals_in_mask) <= 2, 'Bad data mask must have only two pixel values (boolean) - 0 and 1 or ' \
'False and True! The given mask for %s contains the values %s.' % (
self.basename, pixelVals_in_mask)
assert pixelVals_in_mask in [[0, 1], [False, True]], 'Found unsupported pixel values in the given bad ' \
'data mask for %s: %s Only the values True, False, 0 ' \
'and 1 are supported. ' % (self.basename, pixelVals_in_mask)
self._arr = ndarray.astype(np.bool)
def _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=None, fillVal=0):
"""
NOTE: asserts that mapBounds have the same projection like the coordinates in arr_gt
......
......@@ -4,20 +4,22 @@ __author__ = "Daniel Scheffler"
import numpy as np
from pandas import DataFrame
from osgeo import gdal_array
try:
from osgeo import gdal
except ImportError:
import gdal # FIXME this will import this __module__
from ...dtypes.conversion import dTypeDic_NumPy2GDALcompatible
from ...compatibility.gdalnumeric import get_gdalnumeric_func
def get_GDAL_ds_inmem(array, gt=None, prj=None, nodata=None):
"""
"""Convert a numpy array into a GDAL dataset.
NOTE: Possibly the data type has to be automatically changed in order ensure GDAL compatibility!
:param array: <numpy.ndarray> in the shape (rows, columns, bands)
:param array: <numpy.ndarray> in the shape (rows, columns, bands)
:param gt:
:param prj:
:param nodata: nodata value to be set
......@@ -28,6 +30,10 @@ def get_GDAL_ds_inmem(array, gt=None, prj=None, nodata=None):
if len(array.shape) == 3:
array = np.rollaxis(array, 2) # rows,cols,bands => bands,rows,cols
# convert data type to GDAL compatible data type
if gdal_array.NumericTypeCodeToGDALTypeCode(array.dtype) is None:
array = array.astype(dTypeDic_NumPy2GDALcompatible[str(np.dtype(array.dtype))])
OpenNumPyArray = get_gdalnumeric_func('OpenNumPyArray')
ds = OpenNumPyArray(array)
if ds is None:
......
Supports Markdown
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