Commit 222d8f50 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

some bugfixes and improvements

geo.raster.reproject:
- moved get_GDAL_ds_inmem() to new module 'gdal'

geo.vector.geometry:
- boxObj.boxMapYX:  added docstring

geo.vector.topology:
- get_footprint_polygon(): added suppression of shapely self intersection warnings (not yet working) and implemented new keyword 'fix_invalid'

geo.coord_calc:
- calc_FullDataset_corner_positions(): added suppression of shapely self intersection warnings (not yet working)

io.raster.GeoArray:
GeoArray:
    - added property ndim
    - revised __getitem__()
    - from_path(): added docstring and support for negative slices
    - added method save()
    - added method dump()
    - show(): added xlim/ylim parameters for viewing array subsets
- added get_GeoArray_from_GDAL_ds()

new module:
- io.raster.gdal
parent 2c88efe8
......@@ -2,6 +2,7 @@
__author__ = "Daniel Scheffler"
import warnings
import numpy as np
# custom
......@@ -57,11 +58,12 @@ def calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=True, algor
xvalues = np.tile(np.arange(mask_1bit.shape[1]), (1, 2))
outline = np.vstack([topbottom, xvalues])[:, mask].T
poly = Polygon(outline).convex_hull
poly = poly.simplify(20, preserve_topology=True)
poly = Polygon(list(set(poly.exterior.coords))).convex_hull # eliminiate duplicates except of 1 (because poly is closed)
unique_corn_YX = sorted(set(poly.exterior.coords), key=lambda x: poly.exterior.coords[:].index(x)) # [(rows,cols),rows,cols]
with warnings.catch_warnings():
warnings.simplefilter('ignore') # FIXME not working
poly = Polygon(outline).convex_hull
poly = poly.simplify(20, preserve_topology=True)
poly = Polygon(list(set(poly.exterior.coords))).convex_hull # eliminiate duplicates except of 1 (because poly is closed)
unique_corn_YX = sorted(set(poly.exterior.coords), key=lambda x: poly.exterior.coords[:].index(x)) # [(rows,cols),rows,cols]
if assert_four_corners or len(unique_corn_YX)==4:
# sort calculated corners like this: [UL, UR, LL, LR]
......@@ -94,8 +96,8 @@ def calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=True, algor
(cols - list(reversed(mask_1bit[last_dataRow, :])).index(True) + 1)]
if True in [abs(np.diff(i)[0]) > 10 for i in
[StartStopRows_in_first_dataCol, StartStopRows_in_last_dataCol, StartStopCols_in_first_dataRow,
StartStopCols_in_last_dataRow]]:
[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 trueDataCornerPos outside of the image.'''
line_N = ((StartStopCols_in_first_dataRow[1], first_dataRow),
(last_dataCol, StartStopRows_in_last_dataCol[0]))
......
......@@ -17,9 +17,10 @@ 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 ..projection import WKT2EPSG, isProjectedOrGeographic
from ..coord_trafo import pixelToLatLon
from ... import GeoArray
from ..projection import WKT2EPSG, isProjectedOrGeographic
from ..coord_trafo import pixelToLatLon
from ... import GeoArray
from ...io.raster.gdal import get_GDAL_ds_inmem
# dictionary to translate GDAL data types (strings) in corresponding numpy data types
......@@ -212,22 +213,6 @@ def warp_ndarray_OLD(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=N
return out_arr, out_gt, out_prj
def get_GDAL_ds_inmem(array, gt, prj):
if len(array.shape) == 3:
array = np.rollaxis(array, 2) # rows,cols,bands => bands,rows,cols
ds = gdalnumeric.OpenNumPyArray(array)
ds.SetGeoTransform(gt)
ds.SetProjection(prj)
return ds
def get_GeoArray_from_GDAL_ds(ds):
arr = gdalnumeric.DatasetReadAsArray(ds)
if len(arr.shape) == 3:
arr = np.swapaxes(np.swapaxes(arr, 0, 2), 0, 1)
return GeoArray(arr, ds.GetGeoTransform(), ds.GetProjection())
def warp_GeoArray(geoArr, **kwargs):
ndarray = geoArr[:]
return GeoArray(*warp_ndarray(ndarray, geoArr.geotransform, geoArr.projection, **kwargs))
......
......@@ -75,6 +75,10 @@ class boxObj(object):
@property
def boxMapYX(self):
"""Returns a list of YX coordinate tuples for all corners in the order UL_YX, UR_YX, LR_YX, LL_YX.
:return: UL_YX, UR_YX, LR_YX, LL_YX
"""
return shapelyBox2BoxYX(self.mapPoly, coord_type='map')
@boxMapYX.setter
......
......@@ -3,6 +3,7 @@ __author__ = "Daniel Scheffler"
import math
import warnings
from shapely.geometry import shape, Polygon, box
from ..coord_trafo import mapYX2imYX
......@@ -28,14 +29,26 @@ def get_overlap_polygon(poly1, poly2, v=False):
return {'overlap poly':None, 'overlap percentage':0, 'overlap area':0}
def get_footprint_polygon(CornerLonLat):
def get_footprint_polygon(CornerLonLat, fix_invalid=False):
""" Converts a list of coordinates into a shapely polygon object.
:param CornerLonLat: a list of coordinate tuples like [[lon,lat], [lon. lat], ..]
in clockwise or counter clockwise order
:return: a shapyly.Polygon() object
:param fix_invalid: fix invalid output polygon by returning its convex hull (somtimes this can be different)
:return: a shapely.Polygon() object
"""
outpoly = Polygon(CornerLonLat)
assert outpoly.is_valid, 'The given coordinates result in an invalid polygon. Check coordinate order.'
if fix_invalid:
with warnings.catch_warnings():
warnings.simplefilter('ignore') # FIXME not working
outpoly = Polygon(CornerLonLat)
if not outpoly.is_valid:
outpoly = outpoly.convex_hull
else:
outpoly = Polygon(CornerLonLat)
assert outpoly.is_valid, 'The given coordinates result in an invalid polygon. Check coordinate order.' \
'Got coordinates %s.' %CornerLonLat
return outpoly
......
......@@ -10,6 +10,8 @@ from matplotlib import pyplot as plt
# custom
from shapely.geometry import Polygon, box
from osgeo import gdal_array
import dill
try:
from osgeo import gdal
from osgeo import gdalnumeric
......@@ -24,6 +26,7 @@ from ...geo.coord_trafo import mapXY2imXY, imXY2mapXY, transform_any_prj
from ...geo.projection import prj_equal
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
class GeoArray(object):
......@@ -35,8 +38,12 @@ class GeoArray(object):
:param geotransform: GDAL geotransform of the given array or file on disk
:param projection: projection of the given array or file on disk as WKT string
(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']
:param bandnames: names of the bands within the input array, e.g. ['mask_1bit', 'mask_clouds'],
(default: ['B1', 'B2', 'B3', ...])
"""
# FIXME implement compatibility to GDAL VRTs
if type(path_or_array) not in [str, np.ndarray, type(self)]:
raise ValueError("%s parameter 'arg' takes only string "
"or np.ndarray types. Got %s." %(self.__class__.__name__,type(path_or_array)))
......@@ -96,6 +103,11 @@ class GeoArray(object):
return self._shape
@property
def ndim(self):
return len(self.shape)
@property
def rows(self):
return self.shape[0]
......@@ -180,32 +192,44 @@ class GeoArray(object):
def __getitem__(self, given):
# TODO check if array cache contains the needed slice and return data from there
# FIXME negative slices do not work
if self.is_inmem:
# this is only the case if GeoArray has been instanced with an array or GeoArray.to_mem() has been executed
if isinstance(given,(int,float)) and self.bands>2:
# return a band
return self.arr[:,:,given]
elif isinstance(given,(tuple,list)) and len(given)==3 and self.bands==1 and isinstance(given[2],(int,float)):
# in case a third dim is requested from 2D-array -> ignore 3rd dim
return self.arr[given[:2]]
if isinstance(given, (int,float,slice)) and self.ndim==3:
# handle 'given' as index for 3rd (bands) dimension
if self.is_inmem:
return self.arr[:, :, given]
else:
return self.arr[given]
else:
if isinstance(given, tuple):
# behave like a numpy array
getitem_params = given
elif isinstance(given, int) or isinstance(given, slice):
# return only one band
getitem_params = [given]
elif isinstance(given, str):
# behave like a dictionary
if self.bandnames:
if given not in self.bandnames: raise ValueError("'%s' is not a known band." % given)
elif isinstance(given, str):
# behave like a dictionary and return the corresponding band
if self.bandnames:
if given not in self.bandnames:
raise ValueError("'%s' is not a known band. Known bands are: %s"
% (given, ', '.join(list(self.bandnames.keys()))))
if self.is_inmem:
return self.arr[:, :, self.bandnames[given]]
else:
getitem_params = [self.bandnames[given]]
else:
raise ValueError
raise ValueError('String indices are only supported if %s has been instanced with bandnames given.'
%self.__class__.__name__)
elif isinstance(given, (tuple, list)) and len(given)==3 and self.ndim==2:
# in case a third dim is requested from 2D-array -> ignore 3rd dim
if self.is_inmem:
return self.arr[given[:2]]
else:
getitem_params = given[:2]
else:
# behave like a numpy array
if self.is_inmem:
return self.arr[given]
else:
getitem_params = [given] if isinstance(given, slice) else given
if not self.is_inmem:
self._arr_cache = self.from_path(self.arg, getitem_params)
return self._arr_cache
......@@ -227,12 +251,13 @@ class GeoArray(object):
# clean array cache in order to avoid cache pickling
self.flush_cache()
return self.__dict__
def __setstate__(self, state):
"""Defines how the attributes of GMS object are unpickled.
NOTE: This method has been implmemnted because otherwise pickled and unplickled instances show recursion errors
NOTE: This method has been implemented because otherwise pickled and unpickled instances show recursion errors
within __getattr__ when requesting any attribute."""
self.__dict__ = state
......@@ -250,14 +275,23 @@ class GeoArray(object):
def from_path(self, path, getitem_params=None):
# type: (str, list) -> np.ndarray
"""Read a GDAL compatible raster image from disk, with respect to the given image position.
:param path: <str> the file path of the image to read
:param getitem_params: <list> a list of slices in the form [row_slice, col_slice, band_slice]
:return out_arr: <np.ndarray> the output array
"""
ds = gdal.Open(path)
if not ds: raise IOError('Error reading image data at %s.' %path)
R, C, B = ds.RasterYSize, ds.RasterXSize, ds.RasterCount
ds = None
# convert getitem_params to subset area to be read
## convert getitem_params to subset area to be read ##
rS, rE, cS, cE, bS, bE, bL = [None] * 7
# populate rS, rE, cS, cE, bS, bE, bL
if getitem_params:
if len(getitem_params) >= 2:
givenR, givenC = getitem_params[:2]
......@@ -281,7 +315,7 @@ class GeoArray(object):
elif isinstance(givenB, int):
bS = givenB
bE = givenB
elif type(givenB) in [tuple, list]:
elif isinstance(givenB,(tuple,list)):
typesInGivenB = [type(i) for i in givenB]
assert len(list(set(typesInGivenB))) == 1, \
'Mixed data types within the list of bands are not supported.'
......@@ -301,8 +335,17 @@ class GeoArray(object):
bE = bE if bE is not None else B - 1
bL = list(range(bS, bE + 1)) if not bL else bL
# 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
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]
# validate subset area bounds to be read
msg = lambda v, idx, sz: '%s is out of bounds for axis %s with size %s' %(v, idx, sz)
msg = lambda v, idx, sz: '%s is out of bounds for axis %s with size %s' %(v, idx, sz) # FIXME numpy raises that error ONLY for the 2nd axis
for val, axIdx, axSize in zip([rS,rE,cS,cE,bS,bE], [0,0,1,1,2,2], [R,R,C,C,B,B]):
if not 0 <= val <= axSize - 1: raise ValueError(msg(val,axIdx,axSize))
......@@ -339,20 +382,34 @@ class GeoArray(object):
return out_arr
def save(self):
raise NotImplementedError
def save(self, out_path, fmt='ENVI', q=False):
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 self.is_inmem:
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
ds = get_GDAL_ds_inmem(out_arr,self.geotransform, self.projection)
gdalnumeric.SaveArray(out_arr, out_path, format=fmt, prototype=ds)
ds = None
else:
src_ds = gdal.Open(self.filePath)
gdal.Translate(out_path, src_ds, format=fmt)
src_ds = None
def dump(self):
raise NotImplementedError
def dump(self, out_path):
with open(out_path,'w') as outF:
dill.dump(self,outF)
def show(self, band=0, figsize=None, interpolation='none', cmap=None):
# TODO implement slice
def show(self, xlim=None, ylim=None, band=0, figsize=None, interpolation='none', cmap=None):
# handle limits
cS, cE = xlim if isinstance(xlim,(tuple,list)) else (0,self.cols-1)
rS, rE = ylim if isinstance(ylim,(tuple,list)) else (0,self.rows-1)
# show image
plt.figure(figsize=figsize)
is_3D = len(self.shape)>2
plt.imshow(self[:,:,band] if is_3D else self[:,:], interpolation=interpolation, cmap=cmap) # FIXME
plt.imshow(self[rS:rE,cS:cE,band], interpolation=interpolation, cmap=cmap) # FIXME
#plt.imshow(self[:30,:30,band] if is_3D else self[:30,:30],interpolation='none')
plt.show()
......@@ -487,6 +544,13 @@ def _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=None, fillVal=0):
return out_arr, out_gt
def get_GeoArray_from_GDAL_ds(ds):
arr = gdalnumeric.DatasetReadAsArray(ds)
if len(arr.shape) == 3:
arr = np.swapaxes(np.swapaxes(arr, 0, 2), 0, 1)
return GeoArray(arr, ds.GetGeoTransform(), ds.GetProjection())
def get_array_at_mapPosOLD(arr, arr_gt, arr_prj, mapBounds, mapBounds_prj, band2get=None, fillVal=0):
# FIXME mapBounds_prj should not be handled as target projection
"""
......
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
import numpy as np
from osgeo import gdalnumeric
def get_GDAL_ds_inmem(array, gt, prj):
if len(array.shape) == 3:
array = np.rollaxis(array, 2) # rows,cols,bands => bands,rows,cols
ds = gdalnumeric.OpenNumPyArray(array)
ds.SetGeoTransform(gt)
ds.SetProjection(prj)
return ds
\ No newline at end of file
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