Commit 12889c79 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

implemented multiprocessing and GCP based warping into warp_ndarray

geo.raster.reproject:
- warp_ndarray:
    - implemented multiprocessing (new keyword CPUs)
    - added progress bar for GDAL translate and GDAL warp
    - added error threshold to GDAL warp
    - added custom options for GDAL warp
    - added quiet mode

io.raster.GeoArray:
- GeoArray:
    - added dummy version of method 'show_map'

processing:
- added new module progress_mon containing functions to view progress bars
parent f9472f19
......@@ -4,6 +4,7 @@ __author__ = "Daniel Scheffler"
import numpy as np
import warnings
import multiprocessing
try:
from osgeo import gdal
from osgeo import gdalnumeric
......@@ -21,6 +22,7 @@ from ..projection import WKT2EPSG, isProjectedOrGeographic
from ..coord_trafo import pixelToLatLon
from ... import GeoArray
from ...io.raster.gdal import get_GDAL_ds_inmem
from ...processing.progress_mon import printProgress
# dictionary to translate GDAL data types (strings) in corresponding numpy data types
......@@ -221,7 +223,7 @@ def warp_GeoArray(geoArr, **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):
out_alpha=False, targetAlignedPixels=False, gcpList=None, options=None, CPUs=1, progress=True, q=False):
"""
:param ndarray:
......@@ -246,10 +248,14 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=(
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),...]
:param options: <str> additional GDAL options as string, e.g. '-nosrcalpha' or '-order'
:param CPUs: <int> number of CPUs to use (default: None, which means 'all CPUs available')
:param progress: <bool> show progress bar (default: True)
:param q: <bool> quiet mode (default: False)
:return:
"""
# TODO test if this function delivers the exact same output like console version, otherwise implment error_threshold=0.125
# 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
......@@ -258,10 +264,12 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=(
'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
CPUs = CPUs if CPUs else multiprocessing.cpu_count()
progressBarTran = (lambda percent01, message, user_data: printProgress(percent01 * 100,
**{'prefix': 'Translating progress', 'suffix': 'Complete', 'barLength': 50})) if progress and not q else None
progressBarWarp = (lambda percent01, message, user_data: printProgress(percent01 * 100,
**{'prefix': 'Warping progress ', 'suffix': 'Complete', 'barLength': 50})) if progress and not q else None
# not yet implemented
cutlineDSName = 'data/cutline.vrt' #'/vsimem/cutline.shp' TODO cutline from OGR datasource. => implement input shapefile or Geopandas dataframe
......@@ -270,37 +278,76 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=(
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'
rpc = [
"HEIGHT_OFF=1466.05894327379",
"HEIGHT_SCALE=144.837606185489",
"LAT_OFF=38.9266809014185",
"LAT_SCALE=-0.108324009570885",
"LINE_DEN_COEFF=1 -0.000392404256440504 -0.0027925489381758 0.000501819414812054 0.00216726134806561 -0.00185617059201599 0.000183834173326118 -0.00290342803717354 -0.00207181007131322 -0.000900223247894285 -0.00132518281680544 0.00165598132063197 0.00681015244696305 0.000547865679631528 0.00516214646283021 0.00795287690785699 -0.000705040639059332 -0.00254360623317078 -0.000291154885056484 0.00070943440010757",
"LINE_NUM_COEFF=-0.000951099635749339 1.41709976082781 -0.939591985038569 -0.00186609235173885 0.00196881101098923 0.00361741523740639 -0.00282449434932066 0.0115361898794214 -0.00276027843825304 9.37913944402154e-05 -0.00160950221565737 0.00754053609977256 0.00461831968713819 0.00274991122620312 0.000689605203796422 -0.0042482778732957 -0.000123966494595151 0.00307976709897974 -0.000563274426174409 0.00049981716767074",
"LINE_OFF=2199.50159296339",
"LINE_SCALE=2195.852519621",
"LONG_OFF=76.0381768085136",
"LONG_SCALE=0.130066683772651",
"SAMP_DEN_COEFF=1 -0.000632078047521022 -0.000544107268758971 0.000172438016778527 -0.00206391734870399 -0.00204445747536872 -0.000715754551621987 -0.00195545265530244 -0.00168532972557267 -0.00114709980708329 -0.00699131177532728 0.0038551339822296 0.00283631282133365 -0.00436885468926666 -0.00381335885955994 0.0018742043611019 -0.0027263909314293 -0.00237054119407013 0.00246374716379501 -0.00121074576302219",
"SAMP_NUM_COEFF=0.00249293151551852 -0.581492592442025 -1.00947448466175 0.00121597346320039 -0.00552825219917498 -0.00194683170765094 -0.00166012459012905 -0.00338315804553888 -0.00152062885009498 -0.000214562164393127 -0.00219914905535387 -0.000662800177832777 -0.00118644828432841 -0.00180061222825708 -0.00364756875260519 -0.00287273485650089 -0.000540077934728493 -0.00166800463003749 0.000201057249109451 -8.49620129025469e-05",
"SAMP_OFF=3300.34602166792",
"SAMP_SCALE=3297.51222987611"
]
WarpMemoryLimit = 0 # not sure if this is the correct keyword??
# get input dataset (in-MEM)
in_ds = get_GDAL_ds_inmem(ndarray, in_gt, in_prj)
# set RPCs
#if rpcList:
# in_ds.SetMetadata(rpc, "RPC")
# transformerOptions = ['RPC_DEM=data/warp_52_dem.tif']
if CPUs is None or CPUs>1:
gdal.SetConfigOption('GDAL_NUM_THREADS', str(CPUs))
# GDAL Translate if needed
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']
)
in_ds = gdal.Translate(
'', in_ds, format='MEM',
outputSRS = get_SRS(out_prj),
GCPs = gcpList,
callback = progressBarTran
)
# NOTE: options = ['SPARSE_OK=YES'] ## => what is that for?
# GDAL Warp
out_ds = gdal.Warp(
'', in_ds, format='MEM',
dstSRS = get_SRS(out_prj),
outputType = get_GDT(out_dtype) if out_dtype else get_GDT(ndarray.dtype),
xRes = out_gsd[0],
yRes = out_gsd[1],
outputBounds = out_bounds,
outputBoundsSRS = get_SRS(out_bounds_prj),
width = out_XYdims[0],
height = out_XYdims[1],
resampleAlg = rspAlg,
srcNodata = in_nodata,
dstNodata = out_nodata,
srcAlpha = in_alpha,
dstAlpha = out_alpha,
options = options if options else [],
#warpOptions = [],
#transformerOptions = [],
targetAlignedPixels = targetAlignedPixels,
tps = True if gcpList else False,
callback = progressBarWarp,
errorThreshold = 0.125, # this is needed to get exactly the same output like the console version of GDAL warp
)
gdal.SetConfigOption('GDAL_NUM_THREADS', None)
if out_ds is None:
raise Exception('Warping Error: ' + gdal.GetLastErrorMsg())
......
......@@ -414,6 +414,14 @@ class GeoArray(object):
plt.show()
def show_map(self):
raise NotImplementedError
from mpl_toolkits.basemap import Basemap
m = Basemap(llcrnrlon=-119, llcrnrlat=22, urcrnrlon=-64, urcrnrlat=49,
projection='lcc', lat_1=33, lat_2=45, lon_0=-95)
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
"""
......
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
import sys
def printProgress (percent, prefix = '', suffix = '', decimals = 1, barLength = 100):
"""
Call in a loop to create terminal progress bar
:param percent: a value between 0 and 100
:param prefix: - Optional : prefix string (Str)
:param suffix: - Optional : suffix string (Str)
:param decimals: - Optional : positive number of decimals in percent complete (Int)
:param barLength: - Optional : character length of bar (Int)
http://stackoverflow.com/questions/3173320/text-progress-bar-in-the-console
"""
formatStr = "{0:." + str(decimals) + "f}"
percents = formatStr.format(percent)
filledLength = int(round(barLength * percent/100))
bar = '█' * filledLength + '-' * (barLength - filledLength)
sys.stdout.write('\r%s |%s| %s%s %s' % (prefix, bar, percents, '%', suffix)),
if percent==100.:
sys.stdout.write('\n')
sys.stdout.flush()
def tqdm_hook(t):
"""
Wraps tqdm instance. Don't forget to close() or __exit__()
the tqdm instance once you're done with it (easiest using `with` syntax).
Example
-------
> with tqdm(...) as t:
... reporthook = my_hook(t)
... urllib.urlretrieve(..., reporthook=reporthook)
http://stackoverflow.com/questions/3173320/text-progress-bar-in-the-console
"""
last_b = [0]
def inner(b=1, bsize=1, tsize=None):
"""
b : int, optional
Number of blocks just transferred [default: 1].
bsize : int, optional
Size of each block (in tqdm units) [default: 1].
tsize : int, optional
Total size (in tqdm units). If [default: None] remains unchanged.
"""
if tsize is not None:
t.total = tsize
t.update((b - last_b[0]) * bsize)
last_b[0] = b
return inner
def printPercentage(i, i_total):
"""Prints a percentage counter from within a loop.
Example:
for i in range(100+1):
time.sleep(0.1)
printPercentage(i)
:param i:
:param i_total:
:return:
http://stackoverflow.com/questions/3173320/text-progress-bar-in-the-console
"""
sys.stdout.write(('='*i)+(''*(i_total-i))+("\r [ %d"%i+"% ] "))
sys.stdout.flush()
\ 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