Commit 903d411c authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

geo.raster.reproject:

- added a new version of warp_ndarray and renamed the old one to warp_ndarray_OLD: much faster than the old version and no issues when warping 3D arrays
- some modifications to warp_ndarray_OLD
- added get_GDAL_ds_inmem()
- added get_GeoArray_from_GDAL_ds()
- added warp_GeoArray(): a function to warp GeoArray objects

geo.vector.coord_trafo:
- added transform_GCPlist()

 geo.vector.projection:
 - WKT2EPSG() now returns an integer instead of a string

io.raster.GeoArray:
- GeoArray:
    - added whitespace assertion to __init__()
    - added properties 'rows', 'cols', 'bands', 'xgsd', 'ygsd', 'box'
    - dtype: bugfix for error when dtype is requested from in_mem-GeoArray
    - __getitem__(): bugfix for returning the wrong array subset when 'given' slice has length 1
    - show(): added interpolation and cmap keywords
    - get_mapPos(): revision
- _clip_array_at_mapPos:
    - Memory errors due to wrong target dimensions are now catched
- get_array_at_mapPos(): major revision

numeric.array:
- added get_outFillZeroSaturated()
parent 16678c74
...@@ -65,7 +65,8 @@ def calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=True, algor ...@@ -65,7 +65,8 @@ def calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=True, algor
if assert_four_corners or len(unique_corn_YX)==4: if assert_four_corners or len(unique_corn_YX)==4:
# sort calculated corners like this: [UL, UR, LL, LR] # sort calculated corners like this: [UL, UR, LL, LR]
assert len(unique_corn_YX)==4, 'Found %s corner coordinates instead of 4. Exiting.' %len(unique_corn_YX) assert len(unique_corn_YX)==4,\
'Found %s unique corner coordinates instead of 4. Exiting.' %len(unique_corn_YX)
get_dist = lambda P1yx ,P2yx : np.sqrt((P1yx[0]-P2yx[0])**2 + (P1yx[1]-P2yx[1])**2) get_dist = lambda P1yx ,P2yx : np.sqrt((P1yx[0]-P2yx[0])**2 + (P1yx[1]-P2yx[1])**2)
getClosest = lambda tgtYX,srcYXs: srcYXs[np.array([get_dist(tgtYX,srcYX) for srcYX in srcYXs]).argmin()] getClosest = lambda tgtYX,srcYXs: srcYXs[np.array([get_dist(tgtYX,srcYX) for srcYX in srcYXs]).argmin()]
...@@ -118,7 +119,7 @@ def calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=True, algor ...@@ -118,7 +119,7 @@ def calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=True, algor
(last_dataRow, dataCol_in_last_dataRow), (dataRow_in_first_dataCol, first_dataCol)] (last_dataRow, dataCol_in_last_dataRow), (dataRow_in_first_dataCol, first_dataCol)]
getClosest = lambda refR, refC: \ getClosest = lambda refR, refC: \
corners[np.array([np.sqrt((refR - c[0]) ** 2 + (refC - c[1]) ** 2) for c in corners]).argmin()] corners[np.array([np.sqrt((refR - c[0]) ** 2 + (refC - c[1]) ** 2) for c in corners]).argmin()] # FIXME this can also result in only 2 corners
out = [getClosest(refR, refC) for refR, refC in [(0,0), (0, cols-1), (rows-1, 0), (rows-1, cols-1)]] out = [getClosest(refR, refC) for refR, refC in [(0,0), (0, cols-1), (rows-1, 0), (rows-1, cols-1)]]
# if UL[0] == 0 and UR[0] == 0: # if UL[0] == 0 and UR[0] == 0:
# envi.save_image(os.path.abspath('./testing/out/__faulty_mask_1bit.hdr'), self.mask_1bit, dtype= 'uint8', interleave='bsq', ext='.bsq',force=True) # envi.save_image(os.path.abspath('./testing/out/__faulty_mask_1bit.hdr'), self.mask_1bit, dtype= 'uint8', interleave='bsq', ext='.bsq',force=True)
......
...@@ -204,4 +204,17 @@ def latLonToPixel(latLonPairs, path_im=None, geotransform=None, projection=None) ...@@ -204,4 +204,17 @@ def latLonToPixel(latLonPairs, path_im=None, geotransform=None, projection=None)
y = (point[0] - gt[3]) / gt[5] y = (point[0] - gt[3]) / gt[5]
# Add the point to our return array # Add the point to our return array
pixelPairs.append([int(x), int(y)]) pixelPairs.append([int(x), int(y)])
return pixelPairs return pixelPairs
\ No newline at end of file
def transform_GCPlist(gcpList, prj_src, prj_tgt):
"""
:param gcpList: <list> list of ground control points in the output coordinate system
to be used for warping, e.g. [gdal.GCP(mapX,mapY,mapZ,column,row),...]
:param prj_src: WKT string
:param prj_tgt: WKT string
:return:
"""
return [gdal.GCP(*transform_any_prj(prj_src, prj_tgt, GCP.GCPX, GCP.GCPY), GCP.GCPZ, GCP.GCPPixel, GCP.GCPLine)
for GCP in gcpList]
\ No newline at end of file
...@@ -18,6 +18,9 @@ def geotransform2mapinfo(gt,prj): ...@@ -18,6 +18,9 @@ def geotransform2mapinfo(gt,prj):
:rtype: list :rtype: list
""" """
if gt is None:
return None
if gt[2]!=0 or gt[4]!=0: # TODO if gt[2]!=0 or gt[4]!=0: # TODO
raise NotImplementedError('Currently rotated datasets are not supported.') raise NotImplementedError('Currently rotated datasets are not supported.')
srs = osr.SpatialReference() srs = osr.SpatialReference()
...@@ -42,6 +45,10 @@ def mapinfo2geotransform(map_info): ...@@ -42,6 +45,10 @@ def mapinfo2geotransform(map_info):
:returns: GDAL GeoTransform, e.g. [249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0] :returns: GDAL GeoTransform, e.g. [249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0]
:rtype: list :rtype: list
""" """
if map_info is None:
return None
# FIXME rotated datasets are currently not supported -> function must return rotation at gt[2] and gt[4] # FIXME rotated datasets are currently not supported -> function must return rotation at gt[2] and gt[4]
if map_info[1]==1 and map_info[2]==1: if map_info[1]==1 and map_info[2]==1:
ULmapX, ULmapY = float(map_info[3]),float(map_info[4]) ULmapX, ULmapY = float(map_info[3]),float(map_info[4])
......
...@@ -107,9 +107,9 @@ def WKT2EPSG(wkt, epsg=os.environ['GDAL_DATA'].replace('/gdal', '/proj/epsg')): ...@@ -107,9 +107,9 @@ def WKT2EPSG(wkt, epsg=os.environ['GDAL_DATA'].replace('/gdal', '/proj/epsg')):
if m: if m:
code = m.group(1) code = m.group(1)
break break
code = code if code else None # match or no match code = int(code) if code else None # match or no match
else: else:
code = ac code = int(ac)
return code return code
......
...@@ -4,6 +4,12 @@ __author__ = "Daniel Scheffler" ...@@ -4,6 +4,12 @@ __author__ = "Daniel Scheffler"
import numpy as np import numpy as np
import warnings import warnings
try:
from osgeo import gdal
from osgeo import gdalnumeric
except ImportError:
import gdal
import gdalnumeric
# custom # custom
import rasterio import rasterio
...@@ -13,9 +19,21 @@ from rasterio.warp import Resampling ...@@ -13,9 +19,21 @@ from rasterio.warp import Resampling
from ..projection import WKT2EPSG, isProjectedOrGeographic from ..projection import WKT2EPSG, isProjectedOrGeographic
from ..coord_trafo import pixelToLatLon from ..coord_trafo import pixelToLatLon
from ... import GeoArray
def warp_ndarray(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None, outUL=None, out_res=None, # 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): out_extent=None, out_dtype=None, rsp_alg=0, in_nodata=None, out_nodata=None, outExtent_within=True):
"""Reproject / warp a numpy array with given geo information to target coordinate system. """Reproject / warp a numpy array with given geo information to target coordinate system.
...@@ -42,6 +60,7 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None, ...@@ -42,6 +60,7 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None,
:return out_gt: warped gdal GeoTransform :return out_gt: warped gdal GeoTransform
:return out_prj: warped projection as WKT string :return out_prj: warped projection as WKT string
""" """
if not ndarray.flags['OWNDATA']: if not ndarray.flags['OWNDATA']:
temp = np.empty_like(ndarray) temp = np.empty_like(ndarray)
temp[:] = ndarray temp[:] = ndarray
...@@ -55,8 +74,10 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None, ...@@ -55,8 +74,10 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None,
inEPSG, outEPSG = [WKT2EPSG(prj) for prj in [in_prj, out_prj]] inEPSG, outEPSG = [WKT2EPSG(prj) for prj in [in_prj, out_prj]]
assert inEPSG, 'Could not derive input EPSG code.' assert inEPSG, 'Could not derive input EPSG code.'
assert outEPSG, 'Could not derive output EPSG code.' assert outEPSG, 'Could not derive output EPSG code.'
assert in_nodata is None or type(in_nodata) in [int, float] assert in_nodata is None or isinstance(in_nodata, (int, float)),\
assert out_nodata is None or type(out_nodata) in [int, float] 'Received invalid input nodata value: %s; type: %s.' % (in_nodata, type(in_nodata))
assert out_nodata is None or isinstance(out_nodata,(int, float)),\
'Received invalid output nodata value: %s; type: %s.' % (out_nodata, type(out_nodata))
src_crs = {'init': 'EPSG:%s' % inEPSG} src_crs = {'init': 'EPSG:%s' % inEPSG}
dst_crs = {'init': 'EPSG:%s' % outEPSG} dst_crs = {'init': 'EPSG:%s' % outEPSG}
...@@ -75,9 +96,8 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None, ...@@ -75,9 +96,8 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None,
out_dtype = np.int16 if str(out_dtype) == 'int8' else out_dtype out_dtype = np.int16 if str(out_dtype) == 'int8' else out_dtype
ndarray = ndarray.astype(np.int16) if str(in_dtype) == 'int8' else ndarray ndarray = ndarray.astype(np.int16) if str(in_dtype) == 'int8' else ndarray
gt2bounds = lambda gt, r, c: [gt[0], gt[3] + r * gt[5], gt[0] + c * gt[1], gt[3]] # left, bottom, right, top
# get dst_transform # get dst_transform
gt2bounds = lambda gt, r, c: [gt[0], gt[3] + r * gt[5], gt[0] + c * gt[1], gt[3]] # left, bottom, right, top
if out_gt is None and out_extent is None: if out_gt is None and out_extent is None:
if outRowsCols: if outRowsCols:
outUL = [in_gt[0], in_gt[3]] if outUL is None else outUL outUL = [in_gt[0], in_gt[3]] if outUL is None else outUL
...@@ -111,58 +131,77 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None, ...@@ -111,58 +131,77 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None,
out_res = tuple(reversed(abs(px_size_LatLon))) out_res = tuple(reversed(abs(px_size_LatLon)))
print('OUT_RES NOCHMAL CHECKEN: ', out_res) print('OUT_RES NOCHMAL CHECKEN: ', out_res)
dst_transform, out_cols, out_rows = rio_calc_transform( dict_rspInt_rspAlg = \
src_crs, dst_crs, cols, rows, left, bottom, right, top, resolution=out_res) # TODO keyword densify_pts=21 does not exist anymore (moved to transform_bounds()) -> could that be a problem? check warp outputs {0: Resampling.nearest, 1: Resampling.bilinear, 2: Resampling.cubic,
3: Resampling.cubic_spline, 4: Resampling.lanczos, 5: Resampling.average, 6: Resampling.mode}
# check if calculated output dimensions correspond to expected ones and correct them if neccessary var1=True
rows_expected = int(round(abs(top - bottom) / out_res[1], 0)) if var1:
cols_expected = int(round(abs(right - left) / out_res[0], 0)) src_transform = rasterio.transform.from_origin(in_gt[0], in_gt[3], in_gt[1], abs(in_gt[5]))
diff_rows_exp_real, diff_cols_exp_real = abs(out_rows - rows_expected), abs(out_cols - cols_expected) print('calc_trafo_args')
if diff_rows_exp_real > 0.1 or diff_cols_exp_real > 0.1: [print(i, '\n') for i in [src_crs, dst_crs, cols, rows, left, bottom, right, top, out_res]]
assert diff_rows_exp_real < 1.1 and diff_cols_exp_real < 1.1, 'warp_ndarray: The output image size ' \ from ...io.raster.GeoArray import GeoArray
'calculated by rasterio is too far away from the expected output image size.' left, right, bottom, top = GeoArray(ndarray, in_gt, in_prj).box.boundsMap
out_rows, out_cols = rows_expected, cols_expected
# fixes an issue where rio_calc_transform() does not return quadratic output image although input parameters
# correspond to a quadratic image and inEPSG equals outEPSG
aff = list(dst_transform) dst_transform, out_cols, out_rows = rio_calc_transform(
out_gt = out_gt if out_gt else (aff[2], aff[0], aff[1], aff[5], aff[3], aff[4]) src_crs, dst_crs, cols, rows, left, bottom, right, top, resolution=out_res)
src_transform = rasterio.transform.from_origin(in_gt[0], in_gt[3], in_gt[1], in_gt[5])
dict_rspInt_rspAlg = \ out_arr = np.zeros((bands, out_rows, out_cols), out_dtype) \
{0: Resampling.nearest, 1: Resampling.bilinear, 2: Resampling.cubic, if len(ndarray.shape) == 3 else np.zeros((out_rows, out_cols), out_dtype)
3: Resampling.cubic_spline, 4: Resampling.lanczos, 5: Resampling.average, 6: Resampling.mode} print(out_res)
[print(i,'\n') for i in [src_transform, src_crs, dst_transform, dst_crs]]
out_arr = np.zeros((bands, out_rows, out_cols), out_dtype) \ rio_reproject(ndarray, out_arr, src_transform=src_transform, src_crs=src_crs, dst_transform=dst_transform,
if len(ndarray.shape) == 3 else np.zeros((out_rows, out_cols), out_dtype) dst_crs=dst_crs, resampling=dict_rspInt_rspAlg[rsp_alg], src_nodata=in_nodata, dst_nodata=out_nodata)
# FIXME direct passing of src_transform and dst_transform results in a wrong output image. Maybe a rasterio-bug? aff = list(dst_transform)
# rio_reproject(ndarray, out_arr, src_transform=src_transform, src_crs=src_crs, dst_transform=dst_transform, out_gt = out_gt if out_gt else (aff[2], aff[0], aff[1], aff[5], aff[3], aff[4])
# dst_crs=dst_crs, resampling=dict_rspInt_rspAlg[rsp_alg]) # FIXME sometimes output dimensions are not exactly as requested (1px difference)
# FIXME indirect passing causes Future warning else:
with warnings.catch_warnings(): dst_transform, out_cols, out_rows = rio_calc_transform(
warnings.simplefilter('ignore') # FIXME supresses: FutureWarning: GDAL-style transforms are deprecated and will not be supported in Rasterio 1.0. src_crs, dst_crs, cols, rows, left, bottom, right, top, resolution=out_res)
try:
#print('INPUTS') # check if calculated output dimensions correspond to expected ones and correct them if neccessary
#print(ndarray.shape, ndarray.dtype, out_arr.shape, out_arr.dtype) rows_expected = int(round(abs(top - bottom) / out_res[1], 0))
#print(in_gt) cols_expected = int(round(abs(right - left) / out_res[0], 0))
#print(src_crs)
#print(out_gt) #diff_rows_exp_real, diff_cols_exp_real = abs(out_rows - rows_expected), abs(out_cols - cols_expected)
#print(dst_crs) #if diff_rows_exp_real > 0.1 or diff_cols_exp_real > 0.1:
#print(dict_rspInt_rspAlg[rsp_alg]) #assert diff_rows_exp_real < 1.1 and diff_cols_exp_real < 1.1, 'warp_ndarray: ' \
#print(in_nodata) # 'The output image size calculated by rasterio is too far away from the expected output image size.'
#print(out_nodata) # out_rows, out_cols = rows_expected, cols_expected
rio_reproject(ndarray, out_arr, # fixes an issue where rio_calc_transform() does not return quadratic output image although input parameters
src_transform=in_gt, src_crs=src_crs, dst_transform=out_gt, dst_crs=dst_crs, # correspond to a quadratic image and inEPSG equals outEPSG
resampling=dict_rspInt_rspAlg[rsp_alg], src_nodata=in_nodata, dst_nodata=out_nodata)
from matplotlib import pyplot as plt aff = list(dst_transform)
#print(out_arr.shape) out_gt = out_gt if out_gt else (aff[2], aff[0], aff[1], aff[5], aff[3], aff[4])
#plt.figure()
#plt.imshow(out_arr[:,:,1]) out_arr = np.zeros((bands, out_rows, out_cols), out_dtype) \
except KeyError: if len(ndarray.shape) == 3 else np.zeros((out_rows, out_cols), out_dtype)
print(in_dtype, str(in_dtype))
print(ndarray.dtype) with warnings.catch_warnings():
warnings.simplefilter('ignore') # FIXME supresses: FutureWarning: GDAL-style transforms are deprecated and will not be supported in Rasterio 1.0.
try:
#print('INPUTS')
#print(ndarray.shape, ndarray.dtype, out_arr.shape, out_arr.dtype)
#print(in_gt)
#print(src_crs)
#print(out_gt)
#print(dst_crs)
#print(dict_rspInt_rspAlg[rsp_alg])
#print(in_nodata)
#print(out_nodata)
[print(i, '\n') for i in [in_gt, src_crs, out_gt, dst_crs]]
rio_reproject(ndarray, out_arr,
src_transform=in_gt, src_crs=src_crs, dst_transform=out_gt, dst_crs=dst_crs,
resampling=dict_rspInt_rspAlg[rsp_alg], src_nodata=in_nodata, dst_nodata=out_nodata)
from matplotlib import pyplot as plt
#print(out_arr.shape)
#plt.figure()
#plt.imshow(out_arr[:,:,1])
except KeyError:
print(in_dtype, str(in_dtype))
print(ndarray.dtype)
# convert output array axis order to GMS axis order [rows,cols,bands] # convert output array axis order to GMS axis order [rows,cols,bands]
out_arr = out_arr if len(ndarray.shape) == 2 else np.swapaxes(np.swapaxes(out_arr, 0, 1), 1, 2) out_arr = out_arr if len(ndarray.shape) == 2 else np.swapaxes(np.swapaxes(out_arr, 0, 1), 1, 2)
...@@ -171,3 +210,125 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None, ...@@ -171,3 +210,125 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None,
out_arr = out_arr[:outRowsCols[0], :outRowsCols[1]] out_arr = out_arr[:outRowsCols[0], :outRowsCols[1]]
return out_arr, out_gt, out_prj 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))
def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=(None, None),
out_bounds=None, out_bounds_prj=None, out_XYdims = (None,None),
rspAlg='near', in_nodata=None, out_nodata=None, in_alpha=False,
out_alpha=False, targetAlignedPixels=False, gcpList=None):
"""
:param ndarray:
:param in_gt:
:param in_prj:
:param out_prj:
:param out_dtype: gdal.DataType
:param out_gsd:
:param out_bounds: [xmin,ymin,xmax,ymax] set georeferenced extents of output file to be created,
e.g. [440720, 3750120, 441920, 3751320])
(in target SRS by default, or in the SRS specified with -te_srs)
:param out_bounds_prj:
:param out_XYdims:
:param rspAlg: <str> Resampling method to use. Available methods are:
near, bilinear, cubic, cubicspline, lanczos, average, mode, max, min, med, q1, q2
:param in_nodata:
:param out_nodata:
:param in_alpha: <bool> Force the last band of a source image to be considered as a source alpha band.
:param out_alpha: <bool> Create an output alpha band to identify nodata (unset/transparent) pixels
:param targetAlignedPixels: (GDAL >= 1.8.0) (target aligned pixels) align the coordinates of the extent
of the output file to the values of the -tr, such that the aligned extent
includes the minimum extent.
:param gcpList: <list> list of ground control points in the output coordinate system
to be used for warping, e.g. [gdal.GCP(mapX,mapY,mapZ,column,row),...]
:return:
"""
# how to implement: https://svn.osgeo.org/gdal/trunk/autotest/utilities/test_gdalwarp_lib.py
assert str(np.dtype(ndarray.dtype)) in dTypeDic_NumPy2GDAL, "Unknown target datatype '%s'." %ndarray.dtype
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))]
outputType = get_GDT(out_dtype) if out_dtype else get_GDT(ndarray.dtype)
xRes, yRes = out_gsd
width, height = out_XYdims
tps = True if gcpList else False
# not yet implemented
cutlineDSName = 'data/cutline.vrt' #'/vsimem/cutline.shp' TODO cutline from OGR datasource. => implement input shapefile or Geopandas dataframe
cutlineLayer = 'cutline'
cropToCutline = False
cutlineSQL = 'SELECT * FROM cutline'
cutlineWhere = '1 = 1'
warpOptions = [] # ['CUTLINE_ALL_TOUCHED=TRUE'] # this is how to implement extra options
callback = None # mycallback function
callback_data = [0]
transformerOptions = ['SRC_SRS=invalid']
options = '-nosrcalpha'
# get input dataset (in-MEM)
in_ds = get_GDAL_ds_inmem(ndarray, in_gt, in_prj)
if gcpList:
in_ds = gdal.Translate('', in_ds, format='MEM', outputSRS=get_SRS(out_prj), GCPs=gcpList)
gdal.GCP()
# warp
out_ds = gdal.Warp('', in_ds, format='MEM',
dstSRS = get_SRS(out_prj),
outputType = outputType,
xRes = xRes,
yRes = yRes,
outputBounds = out_bounds,
outputBoundsSRS = get_SRS(out_bounds_prj),
width = width,
height = height,
resampleAlg = rspAlg,
srcNodata = in_nodata,
dstNodata = out_nodata,
srcAlpha = in_alpha,
dstAlpha = out_alpha,
warpOptions = warpOptions,
targetAlignedPixels = targetAlignedPixels,
tps = tps,
#transformerOptions=['order=0']
)
if out_ds is None:
raise Exception('Warping Error: ' + gdal.GetLastErrorMsg())
# get outputs
out_arr = gdalnumeric.DatasetReadAsArray(out_ds)
if len(out_arr.shape) == 3:
out_arr = np.swapaxes(np.swapaxes(out_arr, 0, 2), 0, 1)
out_gt = out_ds.GetGeoTransform()
out_prj = out_ds.GetProjection()
# cleanup
in_ds = out_ds = None
return out_arr, out_gt, out_prj
...@@ -9,7 +9,7 @@ from ..coord_trafo import mapYX2imYX ...@@ -9,7 +9,7 @@ from ..coord_trafo import mapYX2imYX
from ..coord_grid import find_nearest_grid_coord from ..coord_grid import find_nearest_grid_coord
def get_overlap_polygon(poly1, poly2, v=True): def get_overlap_polygon(poly1, poly2, v=False):
""" Returns a dict with the overlap of two shapely.Polygon() objects, the overlap percentage and the overlap area. """ Returns a dict with the overlap of two shapely.Polygon() objects, the overlap percentage and the overlap area.
:param poly1: first shapely.Polygon() object :param poly1: first shapely.Polygon() object
...@@ -35,7 +35,7 @@ def get_footprint_polygon(CornerLonLat): ...@@ -35,7 +35,7 @@ def get_footprint_polygon(CornerLonLat):
:return: a shapyly.Polygon() object :return: a shapyly.Polygon() object
""" """
outpoly = Polygon(CornerLonLat) outpoly = Polygon(CornerLonLat)
assert outpoly. is_valid, 'The given coordinates result in an invalid polygon. Check coordinate order.' assert outpoly.is_valid, 'The given coordinates result in an invalid polygon. Check coordinate order.'
return outpoly return outpoly
...@@ -60,6 +60,7 @@ def get_smallest_shapelyImPolyOnGrid_that_contains_shapelyImPoly(shapelyPoly): ...@@ -60,6 +60,7 @@ def get_smallest_shapelyImPolyOnGrid_that_contains_shapelyImPoly(shapelyPoly):
return box(int(xmin),int(ymin), math.ceil(xmax), math.ceil(ymax)) # round min coords off and max coords on return box(int(xmin),int(ymin), math.ceil(xmax), math.ceil(ymax)) # round min coords off and max coords on
def find_line_intersection_point(line1, line2): def find_line_intersection_point(line1, line2):
xdiff = (line1[0][0] - line1[1][0], line2[0][0] - line2[1][0]) xdiff = (line1[0][0] - line1[1][0], line2[0][0] - line2[1][0])
ydiff = (line1[0][1] - line1[1][1], line2[0][1] - line2[1][1]) ydiff = (line1[0][1] - line1[1][1], line2[0][1] - line2[1][1])
......
...@@ -22,9 +22,8 @@ from ...geo.coord_calc import get_corner_coordinates, calc_FullDataset_cor ...@@ -22,9 +22,8 @@ from ...geo.coord_calc import get_corner_coordinates, calc_FullDataset_cor
from ...geo.coord_grid import snap_bounds_to_pixGrid from ...geo.coord_grid import snap_bounds_to_pixGrid
from ...geo.coord_trafo import mapXY2imXY, imXY2mapXY, transform_any_prj from ...geo.coord_trafo import mapXY2imXY, imXY2mapXY, transform_any_prj
from ...geo.projection import prj_equal from ...geo.projection import prj_equal
from ...geo.vector.topology import get_overlap_polygon from ...geo.vector.topology import get_overlap_polygon, get_footprint_polygon
from ...geo.raster.reproject import warp_ndarray from ...geo.vector.geometry import boxObj
class GeoArray(object): class GeoArray(object):
...@@ -52,6 +51,8 @@ class GeoArray(object): ...@@ -52,6 +51,8 @@ class GeoArray(object):
self._dtype = None self._dtype = None
if isinstance(self.arg, str): 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): if not os.path.exists(self.filePath):
raise FileNotFoundError(self.filePath) raise FileNotFoundError(self.filePath)
...@@ -95,10 +96,27 @@ class GeoArray(object): ...@@ -95,10 +96,27 @@ class GeoArray(object):
return self._shape return self._shape
@property
def rows(self):
return self.shape[0]
@property
def cols(self):
return self.shape[1]
@property
def bands(self):
return self.shape[2] if len(self.shape)>2 else 1
@property @property
def dtype(self): def dtype(self):
if self._dtype: if self._dtype:
return self._dtype return self._dtype
elif self.is_inmem:
return self.arr.dtype
else: else:
self.set_gdalDataset_meta() self.set_gdalDataset_meta()
return self._dtype return self._dtype
...@@ -124,6 +142,16 @@ class GeoArray(object): ...@@ -124,6 +142,16 @@ class GeoArray(object):
self._geotransform = gt self._geotransform = gt
@property
def xgsd(self):
return self.geotransform[1]
@property
def ygsd(self):
return abs(self.geotransform[5])
@property @property
def projection(self): def projection(self):
if self._projection: if self._projection:
...@@ -144,12 +172,25 @@ class GeoArray(object): ...@@ -144,12 +172,25 @@ class GeoArray(object):
self._projection = prj self._projection = prj
@property
def box(self):
mapPoly = get_footprint_polygon(get_corner_coordinates(gt=self.geotransform, cols=self.cols,rows=self.rows))
return boxObj(gt=self.geotransform, prj=self.projection, mapPoly=mapPoly)
def __getitem__(self, given): def __getitem__(self, given):
# TODO check if array cache contains the needed slice and return data from there # TODO check if array cache contains the needed slice and return data from there
# FIXME negative slices do not work # FIXME negative slices do not work
if self.is_inmem: if