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

implemented point-wise read processes in GeoArray

geo,vector.topology:
- fill_holes_within_poly(): bugfix for not checking if polygon has holes in case input poly geom_type=='Polygon'

geo.coord_trafo:
- mapXY2imXY(): now also accepts 2D coordinate pair arrays
- imXY2mapXY(): now also accepts 2D coordinate pair arrays

io.raster.GeoArray.GeoArray:
- added function read_pointData(): allows vectorizes point-wise read processes
- __getitem__(): bugfix for wrong return value in case the elements of 'given' contain lists -> allows vectorizes point-wise read processes
parent 82509a5a
...@@ -130,22 +130,40 @@ def mapXY2imXY(mapXY, gt): ...@@ -130,22 +130,40 @@ def mapXY2imXY(mapXY, gt):
# type: (tuple,list) -> tuple # type: (tuple,list) -> tuple
"""Translates given geo coordinates into pixel locations according to the given image geotransform. """Translates given geo coordinates into pixel locations according to the given image geotransform.
:param mapXY: <tuple> The geo coordinates to be translated in the form (x,y). :param mapXY: <tuple, np.ndarray> The geo coordinates to be translated in the form (x,y) or as np.ndarray [Nx1].
:param gt: <list> GDAL geotransform :param gt: <list> GDAL geotransform
:returns: <tuple> image coordinate tuple X/Y (column, row) :returns: <tuple, np.ndarray> image coordinate tuple X/Y (column, row) or np.ndarray [Nx2]
""" """
return (mapXY[0]-gt[0])/gt[1],(mapXY[1]-gt[3])/gt[5] if isinstance(mapXY, np.ndarray):
if mapXY.ndim==1:
mapXY = mapXY.reshape(1,2)
assert mapXY.shape[1]==2, 'An array in shape [Nx2] is needed. Got shape %s.' %mapXY.shape
imXY = np.empty_like(mapXY)
imXY[:,0] = (mapXY[:,0] - gt[0])/gt[1]
imXY[:,1] = (mapXY[:,1] - gt[3])/gt[5]
return imXY
else:
return (mapXY[0]-gt[0])/gt[1],(mapXY[1]-gt[3])/gt[5]
def imXY2mapXY(imXY, gt): def imXY2mapXY(imXY, gt):
# type: (tuple,list) -> tuple # type: (tuple,list) -> tuple
"""Translates given pixel locations into geo coordinates according to the given image geotransform. """Translates given pixel locations into geo coordinates according to the given image geotransform.
:param imXY: <tuple> The image coordinates to be translated in the form (x,y). :param imXY: <tuple, np.ndarray> The image coordinates to be translated in the form (x,y) or as np.ndarray [Nx1].
:param gt: <list> GDAL geotransform :param gt: <list> GDAL geotransform
:returns: <tuple> geo coordinate tuple X/Y (mapX, mapY) :returns: <tuple, np.ndarray> geo coordinate tuple X/Y (mapX, mapY) or np.ndarray [Nx2]
""" """
return (gt[0] + imXY[0] * abs(gt[1])), (gt[3] - imXY[1] * abs(gt[5])) if isinstance(imXY, np.ndarray):
if imXY.ndim==1:
imXY = imXY.reshape(1,2)
assert imXY.shape[1]==2, 'An array in shape [Nx2] is needed. Got shape %s.' %imXY.shape
imXY = np.empty_like(imXY)
imXY[:,0] = gt[0] + imXY[:,0] * abs(gt[1])
imXY[:,1] = gt[3] - imXY[:,1] * abs(gt[5])
return imXY
else:
return (gt[0] + imXY[0] * abs(gt[1])), (gt[3] - imXY[1] * abs(gt[5]))
mapYX2imYX = lambda mapYX,gt: ((mapYX[0]-gt[3])/gt[5], (mapYX[1]-gt[0])/gt[1]) mapYX2imYX = lambda mapYX,gt: ((mapYX[0]-gt[3])/gt[5], (mapYX[1]-gt[0])/gt[1])
......
...@@ -116,7 +116,8 @@ def fill_holes_within_poly(poly): ...@@ -116,7 +116,8 @@ def fill_holes_within_poly(poly):
:return: :return:
""" """
if poly.geom_type == 'Polygon': if poly.geom_type == 'Polygon':
return poly # return only the exterior polygon
return Polygon(poly.exterior)
elif poly.geom_type == 'MultiPolygon': elif poly.geom_type == 'MultiPolygon':
gdf = GeoDataFrame(columns=['geometry']) gdf = GeoDataFrame(columns=['geometry'])
...@@ -125,6 +126,7 @@ def fill_holes_within_poly(poly): ...@@ -125,6 +126,7 @@ def fill_holes_within_poly(poly):
# get the area of each polygon of the multipolygon EXCLUDING the gaps in it # get the area of each polygon of the multipolygon EXCLUDING the gaps in it
gdf['area_filled'] = gdf.apply( gdf['area_filled'] = gdf.apply(
lambda GDF_row: Polygon(np.swapaxes(np.array(GDF_row.geometry.exterior.coords.xy), 0, 1)).area, axis=1) lambda GDF_row: Polygon(np.swapaxes(np.array(GDF_row.geometry.exterior.coords.xy), 0, 1)).area, axis=1)
largest_poly_filled = gdf.ix[gdf['area_filled'].idxmax()]['geometry'] largest_poly_filled = gdf.ix[gdf['area_filled'].idxmax()]['geometry']
# return the outer boundary of the largest polygon # return the outer boundary of the largest polygon
......
...@@ -376,6 +376,11 @@ class GeoArray(object): ...@@ -376,6 +376,11 @@ class GeoArray(object):
%self.__class__.__name__) %self.__class__.__name__)
elif isinstance(given, (tuple, list)) and len(given)==3 and self.ndim==2: elif isinstance(given, (tuple, list)) and len(given)==3 and self.ndim==2:
# handle requests like geoArr[[1,2],[3,4] -> not implemented in from_path if array is not in mem
types = [type(i) for i in given]
if list in types or tuple in types:
self.to_mem()
# in case a third dim is requested from 2D-array -> ignore 3rd dim # in case a third dim is requested from 2D-array -> ignore 3rd dim
if self.is_inmem: if self.is_inmem:
return self.arr[given[:2]] return self.arr[given[:2]]
...@@ -383,6 +388,12 @@ class GeoArray(object): ...@@ -383,6 +388,12 @@ class GeoArray(object):
getitem_params = given[:2] getitem_params = given[:2]
else: else:
if isinstance(given, (tuple, list)):
# handle requests like geoArr[[1,2],[3,4] -> not implemented in from_path if array is not in mem
types = [type(i) for i in given]
if list in types or tuple in types:
self.to_mem()
# behave like a numpy array # behave like a numpy array
if self.is_inmem: if self.is_inmem:
return self.arr[given] return self.arr[given]
...@@ -516,6 +527,7 @@ class GeoArray(object): ...@@ -516,6 +527,7 @@ class GeoArray(object):
elif type(givenB) in [str]: elif type(givenB) in [str]:
bL = [self.bandnames[givenB]] bL = [self.bandnames[givenB]]
# set defaults for not given values # set defaults for not given values
rS = rS if rS is not None else 0 rS = rS if rS is not None else 0
rE = rE if rE is not None else R - 1 rE = rE if rE is not None else R - 1
...@@ -532,7 +544,7 @@ class GeoArray(object): ...@@ -532,7 +544,7 @@ class GeoArray(object):
cE = cE if cE >= 0 else self.columns + cE cE = cE if cE >= 0 else self.columns + cE
bS = bS if bS >= 0 else self.bands + bS bS = bS if bS >= 0 else self.bands + bS
bE = bE if bE >= 0 else self.bands + bE bE = bE if bE >= 0 else self.bands + bE
bL = [b if b >= 0 else (self.bands + b) for b in bL] bL = [b if b >= 0 else (self.bands + b) for b in bL]
# validate subset area bounds to be read # 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) # FIXME numpy raises that error ONLY for the 2nd axis 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
...@@ -881,6 +893,35 @@ class GeoArray(object): ...@@ -881,6 +893,35 @@ class GeoArray(object):
return sub_arr, sub_gt, sub_prj return sub_arr, sub_gt, sub_prj
def read_pointData(self, mapXY_points, mapXY_points_prj=None, band=None):
"""Returns the array values for the given set of X/Y coordinates.
NOTE: If GeoArray has been instanced with a file path, the function will read the dataset into memory.
:param mapXY_points: <np.ndarray, tuple> X/Y coordinates of the points of interest. If a numpy array is
given, it must have the shape [Nx2]
:param mapXY_points_prj: <str, int> WKT string or EPSG code of the projection corresponding to the given
coordinates.
:param band: <int> the band index of the band of interest. If None, the values of all bands are
returned.
:return: np.ndarray with shape [Nx2xbands]
"""
mapXY = mapXY_points if isinstance(mapXY_points, np.ndarray) else np.array(mapXY_points).reshape(1,2)
prj = mapXY_points_prj if mapXY_points_prj else self.prj
assert prj, 'A projection is needed for returning image DNs at specific map X/Y coordinates!'
if not prj_equal(prj1=prj, prj2=self.prj):
mapXY = transform_any_prj(prj, self.prj, mapXY[:,0], mapXY[:,1])
imXY = mapXY2imXY(mapXY, self.geotransform)
imYX = np.fliplr(np.array(imXY)).astype(np.int16)
if imYX.size==2: # only one coordinate pair
Y,X = imYX[0].tolist()
return self[Y,X,band]
else: # multiple coordinate pairs
return self[imYX.T.tolist()+[band]]
def to_mem(self): def to_mem(self):
"""Reads the whole dataset into memory and sets self.arr to the read data. """Reads the whole dataset into memory and sets self.arr to the read data.
""" """
......
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