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

added a lot of feature improvements and further developments

geo.raster
- __init__.py: added __all__
- added module conversion containing new function raster2polygon for polygonizing raster arrays

io.raster:
- gdal:
    - get_GDAL_ds_inmem(): gt and prj are now optional keywords
- GeoArray:
    - added function _alias_property()
    - GeoArray:
        - added attribute basename
        - __init__(): added FileNotFoundError for Python 3
        - added 'nodata' property that return nodata value
        - refactored cols attribute to 'columns'
        - added alias attributes 'cols', 'gt', 'prj'
        - added 'mask_nodata' property that returns a no data mask (uint8)
        - added 'footprint_poly' property that returns a shapely polygon with the corresponding footprint
        - added nodata value derival to set_gdalDataset_meta()
        - save(): added automatic creation of output directory
        - _get_plottable_image(): now respects nodata value
        - show(): now respects nodata value
        - show_map(): now respects nodata value
        - added show_footprint() function that shows a zoomable map containing the footprint polygon as overlay (Jupyter Notebook internal)

io.temp_io => refactored to 'pathgen':
- added function get_generic_outpath()

numeric.array:
- added function find_noDataVal() to derive nodata value from a given numpy array

ptds.__init__:
- updated __version__
parent 6e31fc27
......@@ -16,7 +16,7 @@ __all__=['io',
'processing',
'GeoArray']
__version__ = '20161029_01'
__version__ = '20161101_01'
__author__='Daniel Scheffler'
# Validate GDAL version
......
......@@ -13,7 +13,7 @@ except ImportError:
import gdalconst
from ..processing.shell import subcall_with_output
from ..io.temp_io import get_tempfile
from ..io.pathgen import get_tempfile
def Warp(destNameOrDestDS, srcDSOrSrcDSTab, options = '', format = 'GTiff',
......
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
from . import conversion
from . import reproject
__all__=['conversion',
'reproject',]
\ No newline at end of file
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
from shapely.wkb import loads
import numpy as np
try:
import gdal
import ogr
import osr
except ImportError:
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from ...io.raster.gdal import get_GDAL_ds_inmem
def raster2polygon(array_or_GeoArray, gt=None, prj=None, exact=True):
"""Calculates a footprint polygon for the given array or GeoArray.
:param array_or_GeoArray:
:param gt:
:param prj:
:param exact:
:return:
"""
from ... import GeoArray
geoArr = array_or_GeoArray if isinstance(array_or_GeoArray, GeoArray) else GeoArray(array_or_GeoArray, gt, prj)
src_ds = get_GDAL_ds_inmem(geoArr.mask_nodata.astype(np.uint8), geoArr.gt, geoArr.prj)
src_band = src_ds.GetRasterBand(1)
# Create a memory OGR datasource to put results in.
mem_drv = ogr.GetDriverByName('Memory')
mem_ds = mem_drv.CreateDataSource('out')
srs = osr.SpatialReference()
srs.ImportFromWkt(geoArr.prj)
mem_layer = mem_ds.CreateLayer('poly', srs, ogr.wkbPolygon)
fd = ogr.FieldDefn('DN', ogr.OFTInteger)
mem_layer.CreateField(fd)
# run the algorithm.
result = gdal.Polygonize(src_band, src_band.GetMaskBand(), mem_layer, 0, ["8CONNECTED=8"] if exact else [])
if result is None:
raise Exception(gdal.GetLastErrorMsg())
# extract polygon
mem_layer.SetAttributeFilter('DN = 1')
element = mem_layer.GetNextFeature() # asserts that
shplyPoly = loads(element.GetGeometryRef().ExportToWkb())
return shplyPoly
......@@ -2,7 +2,7 @@
__author__ = "Daniel Scheffler"
from . import raster
from . import temp_io
from . import pathgen
__all__=['raster',
'temp_io']
\ No newline at end of file
'pathgen']
\ No newline at end of file
......@@ -14,4 +14,32 @@ def get_tempfile(ext=None,prefix=None,tgt_dir=None):
prefix = 'py_tools_ds__' if prefix is None else prefix
fd, path = tempfile.mkstemp(prefix=prefix,suffix=ext,dir=tgt_dir)
os.close(fd)
return path
\ No newline at end of file
return path
def get_generic_outpath(dir_out='', fName_out='', prefix='', ext='', create_outDir=True,
prevent_overwriting=False):
"""Generates an output path accourding to the given parameters.
:param dir_out: output directory
:param fName_out: output filename
:param prefix: a prefix for the output filename. ignored if fName_out is given
:param ext: the file extension to use
:param create_outDir: whether to automatically create the output directory or not
:param prevent_overwriting: whether to prevent that a output filename is chosen that already exist in the filesystem
:return:
"""
dir_out = dir_out if dir_out else os.path.abspath(os.path.curdir)
if create_outDir and not os.path.isdir(dir_out): os.makedirs(dir_out)
if not fName_out:
fName_out = '%soutput' %prefix + ('.%s' %ext if ext else '')
if prevent_overwriting:
count=1
while os.path.exists(os.path.join(dir_out, fName_out)):
if count==1:
fName_out += str(count)
else:
fName_out[-1] = str(count)
return os.path.join(dir_out, fName_out)
\ No newline at end of file
......@@ -21,18 +21,27 @@ except ImportError:
import gdalnumeric
from ...geo.coord_calc import get_corner_coordinates, calc_FullDataset_corner_positions
from ...geo.coord_grid import snap_bounds_to_pixGrid
from ...geo.coord_trafo import mapXY2imXY, imXY2mapXY, transform_any_prj
from ...geo.projection import prj_equal, WKT2EPSG, EPSG2WKT
from ...geo.vector.topology import get_overlap_polygon, get_footprint_polygon
from ...geo.vector.geometry import boxObj
from ...io.raster.gdal import get_GDAL_ds_inmem
from ...geo.coord_calc import get_corner_coordinates, calc_FullDataset_corner_positions
from ...geo.coord_grid import snap_bounds_to_pixGrid
from ...geo.coord_trafo import mapXY2imXY, imXY2mapXY, transform_any_prj, reproject_shapelyGeometry
from ...geo.projection import prj_equal, WKT2EPSG, EPSG2WKT
from ...geo.raster.conversion import raster2polygon
from ...geo.vector.topology import get_overlap_polygon, get_footprint_polygon
from ...geo.vector.geometry import boxObj
from ...io.raster.gdal import get_GDAL_ds_inmem
from ...numeric.array import find_noDataVal
def _alias_property(key):
return property(
lambda self: getattr(self, key),
lambda self, val: setattr(self, key, val),
lambda self: delattr(self, key))
class GeoArray(object):
# TODO automatic nodataVal detection => add that value to __init__, adjust calls within CoRegSat
def __init__(self, path_or_array, geotransform=None, projection=None, bandnames=None):
def __init__(self, path_or_array, geotransform=None, projection=None, bandnames=None, nodata=None):
# type: (Any, tuple, str, list) -> GeoArray
"""
......@@ -42,6 +51,7 @@ class GeoArray(object):
(only needed if GeoArray is instanced with an array)
:param bandnames: names of the bands within the input array, e.g. ['mask_1bit', 'mask_clouds'],
(default: ['B1', 'B2', 'B3', ...])
:param nodata nodata value
"""
# FIXME implement compatibility to GDAL VRTs
......@@ -50,20 +60,27 @@ class GeoArray(object):
raise ValueError("%s parameter 'arg' takes only string "
"or np.ndarray types. Got %s." %(self.__class__.__name__,type(path_or_array)))
self.arg = path_or_array
self.arr = self.arg if isinstance(self.arg, np.ndarray) else None
self.filePath = self.arg if isinstance(self.arg, str) and self.arg else None
self._arr_cache = None
self._geotransform = None
self._projection = None
self._shape = None
self._dtype = None
self.arg = path_or_array
self.arr = self.arg if isinstance(self.arg, np.ndarray) else None
self.filePath = self.arg if isinstance(self.arg, str) and self.arg else None
self.basename = os.path.splitext(os.path.basename(self.filePath))[0] if not self.is_inmem else 'IN_MEM'
self._arr_cache = None
self._geotransform = None
self._projection = None
self._shape = None
self._dtype = None
self._nodata = nodata
self._mask_nodata = None
self._footprint_poly = None
if isinstance(self.arg, str):
assert ' ' not in self.arg, "The given path contains whitespaces. This is not supported by GDAL."
if not os.path.exists(self.filePath):
raise IOError(self.filePath) # FileNotFoundError is not available in Python 2.7
try:
raise FileNotFoundError(self.filePath)
except NameError: # FileNotFoundError is not available in Python 2.7
raise IOError(self.filePath)
ds = gdal.Open(self.filePath)
if not ds:
......@@ -86,9 +103,9 @@ class GeoArray(object):
self.bandnames = {'B%s' % band: i for i, band in enumerate(range(1, bands + 1))}
if geotransform:
self.geotransform = geotransform
self.geotransform = geotransform # use property in order to validate given value
if projection:
self.projection = projection
self.projection = projection # use property in order to validate given value
@property
......@@ -116,10 +133,13 @@ class GeoArray(object):
@property
def cols(self):
def columns(self):
return self.shape[1]
cols = _alias_property('columns')
@property
def bands(self):
return self.shape[2] if len(self.shape)>2 else 1
......@@ -159,6 +179,9 @@ class GeoArray(object):
self._geotransform = gt
gt = _alias_property('geotransform')
@property
def xgsd(self):
return self.geotransform[1]
......@@ -189,6 +212,9 @@ class GeoArray(object):
self._projection = prj
prj = _alias_property('projection')
@property
def epsg(self):
return WKT2EPSG(self.projection)
......@@ -201,10 +227,54 @@ class GeoArray(object):
@property
def box(self):
mapPoly = get_footprint_polygon(get_corner_coordinates(gt=self.geotransform, cols=self.cols,rows=self.rows))
mapPoly = get_footprint_polygon(get_corner_coordinates(gt=self.geotransform, cols=self.columns, rows=self.rows))
return boxObj(gt=self.geotransform, prj=self.projection, mapPoly=mapPoly)
@property
def nodata(self):
if self._nodata is not None:
return self._nodata
else:
# try to get nodata value from file
if not self.is_inmem:
self.set_gdalDataset_meta()
if self._nodata is None:
self._nodata = find_noDataVal(self)
if self._nodata == 'ambiguous':
warnings.warn('Nodata value could not be clearly identified. It has been set to None.')
self._nodata = None
return self._nodata
@nodata.setter
def nodata(self, value):
self._nodata = value
@property
def mask_nodata(self):
if self._mask_nodata is not None:
return self._mask_nodata
else:
if self.nodata is not None:
self._mask_nodata = np.where(self[:]==self.nodata,0,1).astype(np.uint8) if self.ndim==2 else \
np.all(np.where(self[:]==self.nodata,0,1), axis=2).astype(np.uint8)
return self._mask_nodata
else:
warnings.warn('Calculation of nodata mask failed due to missing no data value.')
return None
@property
def footprint_poly(self):
if self._footprint_poly is not None:
return self._footprint_poly
else:
self._footprint_poly = raster2polygon(self, exact=False)
return self._footprint_poly
def __getitem__(self, given):
# TODO check if array cache contains the needed slice and return data from there
......@@ -286,7 +356,9 @@ class GeoArray(object):
self._dtype = gdal_array.GDALTypeCodeToNumericTypeCode(ds.GetRasterBand(1).DataType)
self._geotransform = ds.GetGeoTransform()
self._projection = ds.GetProjection()
ds = None
band = ds.GetRasterBand(1)
self._nodata = band.GetNoDataValue() # FIXME this does not support different nodata values within the same file
ds = band = None
def from_path(self, path, getitem_params=None):
......@@ -353,8 +425,8 @@ class GeoArray(object):
# convert negative to positive ones
rS = rS if rS >= 0 else self.rows + rS
rE = rE if rE >= 0 else self.rows + rE
cS = cS if cS >= 0 else self.cols + cS
cE = cE if cE >= 0 else self.cols + cE
cS = cS if cS >= 0 else self.columns + cS
cE = cE if cE >= 0 else self.columns + cE
bS = bS if bS >= 0 else self.bands + bS
bE = bE if bE >= 0 else self.bands + bE
bL = [b if b >= 0 else (self.bands + b) for b in bL]
......@@ -401,6 +473,8 @@ class GeoArray(object):
if not q: print('Writing GeoArray of size %s to %s.' %(self.shape, out_path))
assert self.ndim in [2,3], 'Only 2D- or 3D arrays are supported.'
if not os.path.isdir(os.path.dirname(out_path)): os.makedirs(os.path.dirname(out_path))
if self.is_inmem:
ds = get_GDAL_ds_inmem(self.arr,self.geotransform, self.projection) # expects rows,columns,bands
out_arr = self.arr if self.ndim == 2 else np.swapaxes(np.swapaxes(self.arr, 0, 2), 1, 2) # rows, columns, bands => bands, rows, columns
......@@ -420,18 +494,19 @@ class GeoArray(object):
def _get_plottable_image(self, xlim=None, ylim=None, band=0, res_factor=None, nodataVal=None, out_prj=None):
# handle limits
cS, cE = xlim if isinstance(xlim, (tuple, list)) else (0, self.cols - 1)
cS, cE = xlim if isinstance(xlim, (tuple, list)) else (0, self.columns - 1)
rS, rE = ylim if isinstance(ylim, (tuple, list)) else (0, self.rows - 1)
image2plot = self[rS:rE, cS:cE, band]
gt, prj = self.geotransform, self.projection
transOpt = ['SRC_METHOD=NO_GEOTRANSFORM'] if tuple(gt) == (0, 1, 0, 0, 0, -1) else None
xdim, ydim = None, None
nodataVal = nodataVal if nodataVal is not None else self.nodata
if res_factor != 1. and image2plot.shape[0] * image2plot.shape[1] > 1e6: # shape > 1000*1000
# sample image down
xdim, ydim = (self.cols * res_factor, self.rows * res_factor) if res_factor else \
tuple((np.array([self.cols, self.rows]) / (np.array([self.cols, self.rows]).max() / 1000))) # normalize
xdim, ydim = (self.columns * res_factor, self.rows * res_factor) if res_factor else \
tuple((np.array([self.columns, self.rows]) / (np.array([self.columns, self.rows]).max() / 1000))) # normalize
xdim, ydim = int(xdim), int(ydim)
if xdim or ydim or out_prj:
......@@ -466,11 +541,12 @@ class GeoArray(object):
"""
# get image to plot
nodataVal = nodataVal if nodataVal is not None else self.nodata
image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, res_factor, nodataVal)
# set color palette
palette = cmap if cmap else plt.cm.gray
if nodataVal is not None: # do not show nodata # add auto-detection # TODO
palette = cmap if cmap else plt.cm.gray
if nodataVal is not None: # do not show nodata
image2plot = np.ma.masked_equal(image2plot, nodataVal)
vmin, vmax = np.percentile(image2plot.compressed(),2), np.percentile(image2plot.compressed(),98)
palette.set_bad('aqua', 0)
......@@ -481,8 +557,8 @@ class GeoArray(object):
# show image
plt.figure(figsize=figsize)
plt.imshow(image2plot, palette, interpolation=interpolation,extent=(0,self.cols,self.rows,0),
vmin=vmin,vmax=vmax,) # compressed excludes nodata values
plt.imshow(image2plot, palette, interpolation=interpolation, extent=(0, self.columns, self.rows, 0),
vmin=vmin, vmax=vmax, ) # compressed excludes nodata values
plt.show()
......@@ -508,6 +584,7 @@ class GeoArray(object):
assert self.projection, 'A projection is needed for a map visualization. Got %s.' %self.projection
# get image to plot
nodataVal = nodataVal if nodataVal is not None else self.nodata
image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, res_factor, nodataVal, 'epsg:4326')
# calculate corner coordinates of plot
......@@ -526,13 +603,15 @@ class GeoArray(object):
palette = cmap if cmap else plt.cm.gray
if nodataVal is not None: # do not show nodata
image2plot = np.ma.masked_equal(image2plot, nodataVal)
vmin, vmax = np.percentile(image2plot.compressed(), 2), np.percentile(image2plot.compressed(), 98)
palette.set_bad('aqua', 0)
else:
vmin, vmax = np.percentile(image2plot, 2), np.percentile(image2plot, 98)
palette.set_over ('1')
palette.set_under('0')
# add image to map (y-axis must be inverted for basemap)
m.imshow(np.flipud(image2plot), palette, interpolation=interpolation,
vmin=np.percentile(image2plot,2),vmax=np.percentile(image2plot,98))
m.imshow(np.flipud(image2plot), palette, interpolation=interpolation, vmin=vmin, vmax=vmax)
# add coordinate grid lines
parallels = np.arange(-90, 90., 0.25)
......@@ -555,6 +634,7 @@ class GeoArray(object):
# TODO debug this function
# get image to plot
nodataVal = nodataVal if nodataVal is not None else self.nodata
image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, res_factor, nodataVal)
# calculate corner coordinates of plot
......@@ -580,14 +660,17 @@ class GeoArray(object):
# set color palette
palette = cmap if cmap else plt.cm.gray
if nodataVal is not None: # do not show nodata
if nodataVal is not None: # do not show nodata
image2plot = np.ma.masked_equal(image2plot, nodataVal)
vmin, vmax = np.percentile(image2plot.compressed(), 2), np.percentile(image2plot.compressed(), 98)
palette.set_bad('aqua', 0)
else:
vmin, vmax = np.percentile(image2plot, 2), np.percentile(image2plot, 98)
palette.set_over('1')
palette.set_under('0')
m.imshow(np.flipud(image2plot), palette, interpolation=interpolation,
vmin=np.percentile(image2plot, 2), vmax=np.percentile(image2plot, 98))
# add image to map (y-axis must be inverted for basemap)
m.imshow(np.flipud(image2plot), palette, interpolation=interpolation, vmin=vmin, vmax=vmax)
# add coordinate grid lines
parallels = np.arange(-90, 90., 0.25)
......@@ -602,6 +685,27 @@ class GeoArray(object):
plt.show()
def show_footprint(self):
"""This method is intended to be called from Jupyter Notebook and shows a web map containing the calculated
footprint of GeoArray."""
try:
import folium, geojson
except ImportError:
folium, geojson = None, None
if not folium or not geojson:
raise ImportError(
"This method requires the libraries 'folium' and 'geojson'. They can be installed with "
"the shell command 'pip install folium geojson'.")
lonlatPoly = reproject_shapelyGeometry(self.footprint_poly, self.epsg, 4326)
m = folium.Map(location=tuple(np.array(lonlatPoly.centroid.coords.xy).flatten())[::-1])
gjs = geojson.Feature(geometry=lonlatPoly, properties={})
folium.GeoJson(gjs).add_to(m)
return m
def get_mapPos(self, mapBounds, mapBounds_prj, bandslist=None, arr_gt=None, arr_prj=None, fillVal=0, rspAlg='near'): # TODO implement slice for indexing bands
"""
......@@ -651,6 +755,7 @@ class GeoArray(object):
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
......
......@@ -14,7 +14,7 @@ except ImportError:
from ...compatibility.gdalnumeric import get_gdalnumeric_func
def get_GDAL_ds_inmem(array, gt, prj):
def get_GDAL_ds_inmem(array, gt=None, prj=None):
"""
:param array: <numpy.ndarray> in the shape (rows, columns, bands)
......@@ -28,8 +28,8 @@ def get_GDAL_ds_inmem(array, gt, prj):
OpenNumPyArray = get_gdalnumeric_func('OpenNumPyArray')
ds = OpenNumPyArray(array)
ds.SetGeoTransform(gt)
ds.SetProjection(prj)
if gt: ds.SetGeoTransform(gt)
if prj: ds.SetProjection(prj)
ds.FlushCache() # Write to disk.
return ds
......
......@@ -16,4 +16,24 @@ def get_outFillZeroSaturated(dtype):
dict_outFill = {'int8':-128, 'uint8':0 , 'int16':-9999, 'uint16':9999 , 'float32':-9999.}
dict_outZero = {'int8':0 , 'uint8':1 , 'int16':0 , 'uint16':1 , 'float32':0.}
dict_outSaturated = {'int8':127 , 'uint8':256, 'int16':32767, 'uint16':65535, 'float32':65535.}
return dict_outFill[dtype], dict_outZero[dtype], dict_outSaturated[dtype]
\ No newline at end of file
return dict_outFill[dtype], dict_outZero[dtype], dict_outSaturated[dtype]
def find_noDataVal(pathIm_or_GeoArray,bandIdx=0,sz=3):
"""tries to derive no data value from homogenious corner pixels within 3x3 windows (by default)
:param pathIm_or_GeoArray:
:param bandIdx:
:param sz: window size in which corner pixels are analysed
"""
from .. import GeoArray
geoArr = pathIm_or_GeoArray if isinstance(pathIm_or_GeoArray, GeoArray) else GeoArray(pathIm_or_GeoArray)
get_mean_std = lambda corner_subset: {'mean':np.mean(corner_subset), 'std':np.std(corner_subset)}
UL = get_mean_std(geoArr[0:sz,0:sz,bandIdx])
UR = get_mean_std(geoArr[0:sz,-sz:,bandIdx])
LR = get_mean_std(geoArr[-sz:,-sz:,bandIdx])
LL = get_mean_std(geoArr[-sz:,0:sz,bandIdx])
possVals = [i['mean'] for i in [UL,UR,LR,LL] if i['std']==0]
# possVals==[]: all corners are filled with data; np.std(possVals)==0: noDataVal clearly identified
return None if possVals==[] else possVals[0] if np.std(possVals)==0 else 'ambiguous'
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