Commit 3c234cd5 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Merge branch 'support/cleanup' into develop

parents 45d40422 4f40a143
......@@ -16,6 +16,9 @@
[![build status](https://gitext.gfz-potsdam.de/danschef/arosics/badges/master/build.svg)](https://gitext.gfz-potsdam.de/danschef/arosics/commits/master)
[![coverage report](https://gitext.gfz-potsdam.de/danschef/arosics/badges/master/coverage.svg)](http://danschef.gitext.gfz-potsdam.de/arosics/coverage/)
[![pypi_status](https://img.shields.io/pypi/v/arosics.svg)](https://pypi.python.org/pypi/arosics)
[![license](https://img.shields.io/pypi/l/arosics.svg)](https://gitext.gfz-potsdam.de/danschef/arosics/blob/master/LICENSE)
[![python versions](https://img.shields.io/pypi/pyversions/arosics.svg)](https://img.shields.io/pypi/pyversions/arosics.svg)
See also the latest [coverage](http://danschef.gitext.gfz-potsdam.de/arosics/coverage/) and the [nosetests](http://danschef.gitext.gfz-potsdam.de/arosics/nosetests_reports/nosetests.html) HTML report.
......@@ -58,8 +61,12 @@ Using [conda](https://conda.io/docs/), the recommended approach is:
# create virtual environment for arosics, this is optional
conda create -y -q --name arosics python=3
source activate arosics
conda install -y -q -c conda-forge numpy gdal scikit-image matplotlib pyproj rasterio shapely
conda install -y -q -c conda-forge pyfftw basemap pykrige # these libraries are optional
conda install -y -q -c conda-forge numpy gdal scikit-image matplotlib pyproj rasterio shapely geopandas
# optional libraries:
conda install -y -q -c conda-forge basemap pykrige
conda install -y -q -c conda-forge pyfftw # Linux and MacOS
conda install -y -q -c jesserobertson pyfftw # Windows
```
To install AROSICS, use the pip installer:
......
......@@ -29,6 +29,10 @@ Status
:target: http://danschef.gitext.gfz-potsdam.de/arosics/coverage/
.. image:: https://img.shields.io/pypi/v/arosics.svg
:target: https://pypi.python.org/pypi/arosics
.. image:: https://img.shields.io/pypi/l/arosics.svg
:target: https://gitext.gfz-potsdam.de/danschef/arosics/blob/master/LICENSE
.. image:: https://img.shields.io/pypi/pyversions/arosics.svg
:target: https://img.shields.io/pypi/pyversions/arosics.svg
See also the latest coverage_ report and the nosetests_ HTML report.
......@@ -52,8 +56,12 @@ Using conda_, the recommended approach is:
# create virtual environment for arosics, this is optional
conda create -y -q --name arosics python=3
source activate arosics
conda install -y -q -c conda-forge numpy gdal scikit-image matplotlib pyproj rasterio shapely
conda install -y -q -c conda-forge pyfftw basemap pykrige # these libraries are optional
conda install -y -q -c conda-forge numpy gdal scikit-image matplotlib pyproj rasterio shapely geopandas
# optional libraries:
conda install -y -q -c conda-forge basemap pykrige
conda install -y -q -c conda-forge pyfftw # Linux and MacOS
conda install -y -q -c jesserobertson pyfftw # Windows
To install AROSICS, use the pip installer:
......
# -*- coding: utf-8 -*-
import os
import re
import shutil
import subprocess
import time
import warnings
from copy import copy
......@@ -25,21 +22,20 @@ from skimage.exposure import rescale_intensity
# internal modules
from .DeShifter import DESHIFTER, _dict_rspAlg_rsp_Int
from . import geometry as GEO
from . import io as IO
from . import plotting as PLT
from geoarray import GeoArray
from py_tools_ds.convenience.object_oriented import alias_property
from py_tools_ds.geo.coord_calc import corner_coord_to_minmax, get_corner_coordinates
from py_tools_ds.geo.coord_calc import get_corner_coordinates
from py_tools_ds.geo.vector.topology import get_overlap_polygon, get_smallest_boxImYX_that_contains_boxMapYX
from py_tools_ds.geo.projection import prj_equal, get_proj4info
from py_tools_ds.geo.vector.geometry import boxObj, round_shapelyPoly_coords
from py_tools_ds.geo.coord_grid import move_shapelyPoly_to_image_grid
from py_tools_ds.geo.coord_trafo import pixelToMapYX, reproject_shapelyGeometry, mapXY2imXY, imXY2mapXY
from py_tools_ds.geo.coord_trafo import reproject_shapelyGeometry, mapXY2imXY, imXY2mapXY
from py_tools_ds.geo.raster.reproject import warp_ndarray
from py_tools_ds.geo.map_info import geotransform2mapinfo
from py_tools_ds.numeric.vector import find_nearest
from py_tools_ds.similarity.raster import calc_ssim
from py_tools_ds.io.vector.writer import write_shp
__author__ = 'Daniel Scheffler'
......@@ -293,9 +289,9 @@ class COREG(object):
self._get_overlap_properties()
if self.v and self.path_verbose_out:
IO.write_shp(os.path.join(self.path_verbose_out, 'poly_imref.shp'), self.ref.poly, self.ref.prj)
IO.write_shp(os.path.join(self.path_verbose_out, 'poly_im2shift.shp'), self.shift.poly, self.shift.prj)
IO.write_shp(os.path.join(self.path_verbose_out, 'overlap_poly.shp'), self.overlap_poly, self.ref.prj)
write_shp(os.path.join(self.path_verbose_out, 'poly_imref.shp'), self.ref.poly, self.ref.prj)
write_shp(os.path.join(self.path_verbose_out, 'poly_im2shift.shp'), self.shift.poly, self.shift.prj)
write_shp(os.path.join(self.path_verbose_out, 'overlap_poly.shp'), self.overlap_poly, self.ref.prj)
# FIXME: transform_mapPt1_to_mapPt2(im2shift_center_map, ds_imref.GetProjection(), ds_im2shift.GetProjection())
# FIXME später basteln für den fall, dass projektionen nicht gleich sind
......@@ -307,8 +303,8 @@ class COREG(object):
self._get_clip_window_properties() # sets self.matchBox, self.otherBox and much more
if self.v and self.path_verbose_out and self.matchBox.mapPoly and self.success is not False:
IO.write_shp(os.path.join(self.path_verbose_out, 'poly_matchWin.shp'),
self.matchBox.mapPoly, self.matchBox.prj)
write_shp(os.path.join(self.path_verbose_out, 'poly_matchWin.shp'),
self.matchBox.mapPoly, self.matchBox.prj)
self.success = False if self.success is False or not self.matchBox.boxMapYX else None
self._coreg_info = None # private attribute to be filled by self.coreg_info property
......@@ -741,8 +737,8 @@ class COREG(object):
print('Target window size %s not possible due to too small overlap area or window position too close '
'to an image edge. New matching window size: %s.' % (self.win_size_XY, match_win_size_XY))
# IO.write_shp('/misc/hy5/scheffler/Temp/matchMapPoly.shp', matchBox.mapPoly,matchBox.prj)
# IO.write_shp('/misc/hy5/scheffler/Temp/otherMapPoly.shp', otherBox.mapPoly,otherBox.prj)
# write_shp('/misc/hy5/scheffler/Temp/matchMapPoly.shp', matchBox.mapPoly,matchBox.prj)
# write_shp('/misc/hy5/scheffler/Temp/otherMapPoly.shp', otherBox.mapPoly,otherBox.prj)
def _get_image_windows_to_match(self):
"""Reads the matching window and the other window using subset read, and resamples the other window to the
......@@ -1410,150 +1406,3 @@ class COREG(object):
q=self.q)
self.deshift_results = DS.correct_shifts()
return self.deshift_results
def _correct_shifts_OLD(self): # pragma: no cover
if self.success:
if not os.path.exists(os.path.dirname(self.path_out)):
os.makedirs(os.path.dirname(self.path_out))
equal_prj = prj_equal(self.ref.prj, self.shift.prj)
if equal_prj and not self.align_grids and not self.match_gsd and \
self.out_gsd in [None, [self.shift.xgsd, self.shift.ygsd]]:
self._shift_image_by_updating_map_info()
elif equal_prj and self.align_grids: # match_gsd and out_gsd are respected
self._align_coordinate_grids()
else: # match_gsd and out_gsd are respected ### TODO: out_proj implementieren
self._resample_without_grid_aligning()
else:
warnings.warn('No result written because detection of image displacements failed.')
def _shift_image_by_updating_map_info(self): # pragma: no cover
if not self.q:
print('\nWriting output...')
ds_im2shift = gdal.Open(self.shift.path)
if not ds_im2shift.GetDriver().ShortName == 'ENVI': # FIXME laaangsam
if self.CPUs is not None and self.CPUs > 1:
IO.convert_gdal_to_bsq__mp(self.shift.path, self.path_out)
else:
os.system('gdal_translate -of ENVI %s %s' % (self.shift.path, self.path_out))
file2getHdr = self.path_out
else:
shutil.copy(self.shift.path, self.path_out)
file2getHdr = self.shift.path
del ds_im2shift
path_hdr = '%s.hdr' % os.path.splitext(file2getHdr)[0]
path_hdr = '%s.hdr' % file2getHdr if not os.path.exists(path_hdr) else path_hdr
path_hdr = None if not os.path.exists(path_hdr) else path_hdr
assert path_hdr, 'No header file found for %s. Applying shifts failed.' % file2getHdr
new_originY, new_originX = pixelToMapYX([self.x_shift_px, self.y_shift_px],
geotransform=self.shift.gt, projection=self.shift.prj)[0]
map_xshift, map_yshift = new_originX - self.shift.gt[0], new_originY - self.shift.gt[3]
if not self.q:
print('Calculated map shifts (X,Y): %s, %s' % (map_xshift, map_yshift))
with open(path_hdr, 'r') as inF:
content = inF.read()
map_info_line = [i for i in content.split('\n') if re.search('map info [\w+]*=', i, re.I)][0]
map_info_string = re.search('{([\w+\s\S]*)}', map_info_line, re.I).group(1)
map_info = [i.strip() for i in map_info_string.split(',')]
if not self.q:
print('Original map info: %s' % map_info)
new_map_info = map_info
new_map_info[3] = str(float(map_info[3]) + map_xshift)
new_map_info[4] = str(float(map_info[4]) + map_yshift)
if not self.q:
print('Updated map info: %s' % new_map_info)
new_map_info_line = 'map info = { %s }' % ' , '.join(new_map_info)
new_content = content.replace(map_info_line, new_map_info_line)
outpath_hdr = '%s.hdr' % os.path.splitext(self.path_out)[0]
with open(outpath_hdr, 'w') as outF:
outF.write(new_content)
if not self.q:
print('\nCoregistered image written to %s.' % self.path_out)
def _align_coordinate_grids(self): # pragma: no cover
xgsd, ygsd = (self.ref.xgsd, self.ref.ygsd) if self.match_gsd else self.out_gsd if self.out_gsd \
else (self.shift.xgsd, self.shift.ygsd) # self.match_gsd overrides self.out_gsd in __init__
if not self.q:
print('Resampling shift image in order to align coordinate grid of shift image to reference image. ')
def is_alignable(gsd1, gsd2):
"""Check if pixel sizes are divisible"""
return max(gsd1, gsd2) % min(gsd1, gsd2) == 0
if not is_alignable(self.ref.xgsd, xgsd) or not is_alignable(self.ref.ygsd, ygsd) and not self.q:
print("\nWARNING: The coordinate grids of the reference image and the image to be shifted cannot be "
"aligned because their pixel sizes are not exact multiples of each other (ref [X/Y]: "
"%s %s; shift [X/Y]: %s %s). Therefore the pixel size of the reference image is chosen for the "
"resampled output image. If you don´t like that you can use the '-out_gsd' parameter to set an "
"appropriate output pixel size.\n" % (self.ref.xgsd, self.ref.ygsd, xgsd, ygsd))
xgsd, ygsd = self.ref.xgsd, self.ref.ygsd
r_xmin, r_xmax, r_ymin, r_ymax = corner_coord_to_minmax(get_corner_coordinates(self.ref.path))
imref_xgrid = np.arange(r_xmin, r_xmax, self.ref.xgsd)
imref_ygrid = np.arange(r_ymin, r_ymax, self.ref.ygsd)
s_xmin, s_xmax, s_ymin, s_ymax = corner_coord_to_minmax(get_corner_coordinates(self.shift.path))
xmin, xmax = find_nearest(imref_xgrid, s_xmin, 'off', extrapolate=1), find_nearest(imref_xgrid, s_xmax, 'on',
extrapolate=1)
ymin, ymax = find_nearest(imref_ygrid, s_ymin, 'off', extrapolate=1), find_nearest(imref_ygrid, s_ymax, 'on',
extrapolate=1)
new_originY, new_originX = pixelToMapYX([self.x_shift_px, self.y_shift_px],
geotransform=self.shift.gt, projection=self.shift.prj)[0]
map_xshift, map_yshift = new_originX - self.shift.gt[0], new_originY - self.shift.gt[3]
if not self.q:
print('Calculated map shifts (X,Y): %s, %s' % (map_xshift, map_yshift))
gt_shifted = list(self.shift.gt[:])
gt_shifted[0], gt_shifted[3] = new_originX, new_originY
pathVRT = os.path.splitext(os.path.basename(self.shift.path))[0] + '.vrt'
ds_im2shift = gdal.Open(self.shift.path)
dst_ds = gdal.GetDriverByName('VRT').CreateCopy(pathVRT, ds_im2shift, 0)
dst_ds.SetGeoTransform(gt_shifted)
del dst_ds, ds_im2shift
cmd = "gdalwarp -r cubic -tr %s %s -t_srs '%s' -of ENVI %s %s -overwrite -te %s %s %s %s" \
% (xgsd, ygsd, self.ref.prj, pathVRT, self.path_out, xmin, ymin, xmax, ymax)
output = subprocess.check_output(cmd, shell=True)
if output != 1 and os.path.exists(self.path_out):
if not self.q:
print('Coregistered and resampled image written to %s.' % self.path_out)
else:
print(output)
self._handle_error(RuntimeError('Resampling failed.'))
def _resample_without_grid_aligning(self): # pragma: no cover
xgsd, ygsd = (self.ref.xgsd, self.ref.ygsd) if self.match_gsd else self.out_gsd if self.out_gsd \
else (self.shift.xgsd, self.shift.ygsd) # self.match_gsd overrides self.out_gsd in __init__
new_originY, new_originX = pixelToMapYX([self.x_shift_px, self.y_shift_px],
geotransform=self.shift.gt, projection=self.shift.prj)[0]
map_xshift, map_yshift = new_originX - self.shift.gt[0], new_originY - self.shift.gt[3]
if not self.q:
print('Calculated map shifts (X,Y): %s, %s' % (map_xshift, map_yshift))
gt_shifted = list(self.shift.gt[:])
gt_shifted[0], gt_shifted[3] = new_originX, new_originY
pathVRT = os.path.splitext(os.path.basename(self.shift.path))[0] + '.vrt'
ds_im2shift = gdal.Open(self.shift.path)
drv = gdal.GetDriverByName('VRT')
dst_ds = drv.CreateCopy(pathVRT, ds_im2shift, 0)
dst_ds.SetGeoTransform(gt_shifted)
del dst_ds, ds_im2shift
if not self.q:
print('Resampling shift image without aligning of coordinate grids. ')
cmd = "gdalwarp -r %s -tr %s %s -t_srs '%s' -of ENVI %s %s -overwrite" \
% (self.rspAlg_calc, xgsd, ygsd, self.ref.prj, pathVRT, self.path_out)
output = subprocess.check_output(cmd, shell=True)
if output != 1 and os.path.exists(self.path_out):
if not self.q:
print('Coregistered and resampled image written to %s.' % self.path_out)
else:
print(output)
self._handle_error(RuntimeError('Resampling failed.'))
......@@ -20,11 +20,11 @@ from skimage.transform import AffineTransform, PolynomialTransform
# internal modules
from .CoReg import COREG
from . import io as IO
from py_tools_ds.geo.projection import isProjectedOrGeographic, isLocal, get_UTMzone, dict_to_proj4, proj4_to_WKT
from py_tools_ds.io.pathgen import get_generic_outpath
from py_tools_ds.processing.progress_mon import ProgressBar
from py_tools_ds.geo.vector.conversion import points_to_raster
from py_tools_ds.io.vector.writer import write_shp
from geoarray import GeoArray
from .CoReg import GeoArray_CoReg # noqa F401 # flake8 issue
......@@ -611,11 +611,11 @@ class Tie_Point_Grid(object):
GDF2pass = GDF if not skip_nodata else GDF[GDF[skip_nodata_col] != self.outFillVal]
# replace boolean values (cannot be written)
GDF2pass = GDF2pass.replace(False, 0) # replace all booleans where column dtype is not np.bool but np.object
GDF2pass = GDF2pass.replace(True, 1)
for col in GDF2pass.columns:
if GDF2pass[col].dtype == np.bool:
GDF2pass[col] = GDF2pass[col].astype(int)
GDF2pass = GDF2pass.replace(False, 0) # replace all remaining booleans where dtype is not np.bool but np.object
GDF2pass = GDF2pass.replace(True, 1)
path_out = path_out if path_out else \
get_generic_outpath(dir_out=os.path.join(self.dir_out, 'CoRegPoints'),
......@@ -637,7 +637,7 @@ class Tie_Point_Grid(object):
fName_out = "CoRegPoints_grid%s_ws%s.shp" % (self.grid_res, self.COREG_obj.win_size_XY)
path_out = os.path.join(self.dir_out, fName_out)
IO.write_shp(path_out, shapely_points, prj=self.COREG_obj.shift.prj, attrDict=attr_dicts)
write_shp(path_out, shapely_points, prj=self.COREG_obj.shift.prj, attrDict=attr_dicts)
def to_vectorfield(self, path_out=None, fmt=None, mode='md'):
# type: (str) -> GeoArray
......@@ -680,45 +680,6 @@ class Tie_Point_Grid(object):
return out_GA
def _to_Raster_using_KrigingOLD(self, attrName, skip_nodata=1, skip_nodata_col='ABS_SHIFT', outGridRes=None,
path_out=None, tilepos=None): # pragma: no cover
warnings.warn(DeprecationWarning("'to_Raster_using_KrigingOLD' is deprecated. Use to_Raster_using_Kriging "
"instead.")) # TODO delete
GDF = self.CoRegPoints_table
GDF2pass = GDF if not skip_nodata else GDF[GDF[skip_nodata_col] != self.outFillVal]
# subset if tilepos is given
rows, cols = tilepos if tilepos else self.shift.shape
GDF2pass = GDF2pass.loc[(GDF2pass['X_IM'] >= cols[0]) & (GDF2pass['X_IM'] <= cols[1]) &
(GDF2pass['Y_IM'] >= rows[0]) & (GDF2pass['Y_IM'] <= rows[1])]
X_coords, Y_coords, ABS_SHIFT = GDF2pass['X_UTM'], GDF2pass['Y_UTM'], GDF2pass[attrName]
xmin, ymin, xmax, ymax = GDF2pass.total_bounds
grid_res = outGridRes if outGridRes else int(min(xmax - xmin, ymax - ymin) / 250)
grid_x, grid_y = np.arange(xmin, xmax + grid_res, grid_res), np.arange(ymax, ymin - grid_res, -grid_res)
# Reference: P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology,
# (Cambridge University Press, 1997) 272 p.
from pykrige.ok import OrdinaryKriging
OK = OrdinaryKriging(X_coords, Y_coords, ABS_SHIFT, variogram_model='spherical', verbose=False)
zvalues, sigmasq = OK.execute('grid', grid_x, grid_y) # ,backend='C',)
path_out = path_out if path_out else \
get_generic_outpath(dir_out=os.path.join(self.dir_out, 'CoRegPoints'),
fName_out="Kriging__%s__grid%s_ws(%s_%s).tif"
% (attrName, self.grid_res, self.COREG_obj.win_size_XY[0],
self.COREG_obj.win_size_XY[1]))
print('Writing %s ...' % path_out)
# add a half pixel grid points are centered on the output pixels
xmin, ymin, xmax, ymax = xmin - grid_res / 2, ymin - grid_res / 2, xmax + grid_res / 2, ymax + grid_res / 2
IO.write_numpy_to_image(zvalues, path_out, gt=(xmin, grid_res, 0, ymax, 0, -grid_res),
prj=self.COREG_obj.shift.prj)
return zvalues
def to_Raster_using_Kriging(self, attrName, skip_nodata=1, skip_nodata_col='ABS_SHIFT', outGridRes=None,
fName_out=None, tilepos=None, tilesize=500, mp=None):
......@@ -752,17 +713,6 @@ class Tie_Point_Grid(object):
GDF = self.CoRegPoints_table
GDF2pass = GDF if not skip_nodata else GDF[GDF[skip_nodata_col] != self.outFillVal]
# # subset if tilepos is given
# # overlap_factor =
# rows,cols = tilepos if tilepos else self.tgt_shape
# xvals, yvals = np.sort(GDF2pass['X_IM'].values.flat),np.sort(GDF2pass['Y_IM'].values.flat)
# cS,cE = UTL.find_nearest(xvals,cols[0],'off',1), UTL.find_nearest(xvals,cols[1],'on',1)
# rS,rE = UTL.find_nearest(yvals,rows[0],'off',1), UTL.find_nearest(yvals,rows[1],'on',1)
# # GDF2pass = GDF2pass.loc[(GDF2pass['X_IM']>=cols[0])&(GDF2pass['X_IM']<=cols[1])&
# # (GDF2pass['Y_IM']>=rows[0])&(GDF2pass['Y_IM']<=rows[1])]
# GDF2pass = GDF2pass.loc[(GDF2pass['X_IM']>=cS)&(GDF2pass['X_IM']<=cE)&
# (GDF2pass['Y_IM']>=rS)&(GDF2pass['Y_IM']<=rE)]
X_coords, Y_coords, ABS_SHIFT = GDF2pass['X_UTM'], GDF2pass['Y_UTM'], GDF2pass[attrName]
xmin, ymin, xmax, ymax = GDF2pass.total_bounds
......@@ -783,11 +733,12 @@ class Tie_Point_Grid(object):
fName_out = fName_out if fName_out else \
"Kriging__%s__grid%s_ws%s.tif" % (attrName, self.grid_res, self.COREG_obj.win_size_XY)
path_out = get_generic_outpath(dir_out=self.dir_out, fName_out=fName_out)
print('Writing %s ...' % path_out)
# add a half pixel grid points are centered on the output pixels
xmin, ymin, xmax, ymax = xmin - grid_res / 2, ymin - grid_res / 2, xmax + grid_res / 2, ymax + grid_res / 2
IO.write_numpy_to_image(zvalues, path_out, gt=(xmin, grid_res, 0, ymax, 0, -grid_res),
prj=self.COREG_obj.shift.prj)
GeoArray(zvalues,
geotransform=(xmin, grid_res, 0, ymax, 0, -grid_res),
projection=self.COREG_obj.shift.prj).save(path_out)
return zvalues
......
......@@ -12,8 +12,8 @@ from arosics.Tie_Point_Grid import Tie_Point_Grid
__author__ = """Daniel Scheffler"""
__email__ = 'daniel.scheffler@gfz-potsdam.de'
__version__ = '0.6.7'
__versionalias__ = '2017-11-15_01'
__version__ = '0.6.8'
__versionalias__ = '2017-11-16_01'
__all__ = ['COREG',
'COREG_LOCAL',
'DESHIFTER',
......
# -*- coding: utf-8 -*-
import ctypes
import multiprocessing
import os
import time
import numpy as np
import ogr
import osr
try:
import gdal
except ImportError:
from osgeo import gdal
from spectral.io import envi
# internal modules
from .utilities import get_image_tileborders, convertGdalNumpyDataType
from py_tools_ds.geo.map_info import geotransform2mapinfo
from py_tools_ds.geo.projection import EPSG2WKT
from py_tools_ds.dtypes.conversion import get_dtypeStr
def wait_if_used(path_file, lockfile, timeout=100, try_kill=0):
globs = globals()
same_gdalRefs = [k for k, v in globs.items() if
isinstance(globs[k], gdal.Dataset) and globs[k].GetDescription() == path_file]
t0 = time.time()
def update_same_gdalRefs(sRs): return [sR for sR in sRs if sR in globals() and globals()[sR] is not None]
while same_gdalRefs != [] or os.path.exists(lockfile):
if os.path.exists(lockfile):
continue
if time.time() - t0 > timeout:
if try_kill:
for sR in same_gdalRefs:
globals()[sR] = None
print('had to kill %s' % sR)
else:
if os.path.exists(lockfile):
os.remove(lockfile)
raise TimeoutError('The file %s is permanently used by another variable.' % path_file)
same_gdalRefs = update_same_gdalRefs(same_gdalRefs)
def write_envi(arr, outpath, gt=None, prj=None):
if gt or prj:
assert gt and prj, 'gt and prj must be provided together or left out.'
meta = {'map info': geotransform2mapinfo(gt, prj), 'coordinate system string': prj} if gt else None
shape = (arr.shape[0], arr.shape[1], 1) if len(arr.shape) == 3 else arr.shape
out = envi.create_image(outpath, metadata=meta, shape=shape, dtype=arr.dtype, interleave='bsq', ext='.bsq',
force=True) # 8bit for multiple masks in one file
out_mm = out.open_memmap(writable=True)
out_mm[:, :, 0] = arr
def wfa(p, c): # pragma: no cover
try:
with open(p, 'a') as of:
of.write(c)
except Exception:
pass
shared_array = None
def init_SharedArray_in_globals(dims):
rows, cols = dims
global shared_array
shared_array_base = multiprocessing.Array(ctypes.c_double, rows * cols)
shared_array = np.ctypeslib.as_array(shared_array_base.get_obj())
shared_array = shared_array.reshape(rows, cols)
def gdal_read_subset(fPath, pos, bandNr):
(rS, rE), (cS, cE) = pos
ds = gdal.Open(fPath)
data = ds.GetRasterBand(bandNr).ReadAsArray(cS, rS, cE - cS + 1, rE - rS + 1)
del ds
return data
def fill_arr(argDict, def_param=shared_array):
pos = argDict.get('pos')
func = argDict.get('func2call')
args = argDict.get('func_args', [])
kwargs = argDict.get('func_kwargs', {})
(rS, rE), (cS, cE) = pos
shared_array[rS:rE + 1, cS:cE + 1] = func(*args, **kwargs)
def gdal_ReadAsArray_mp(fPath, bandNr, tilesize=1500):
ds = gdal.Open(fPath)
rows, cols = ds.RasterYSize, ds.RasterXSize
del ds
init_SharedArray_in_globals((rows, cols))
tilepos = get_image_tileborders([tilesize, tilesize], (rows, cols))
fill_arr_argDicts = [{'pos': pos, 'func2call': gdal_read_subset, 'func_args': (fPath, pos, bandNr)} for pos in
tilepos]
with multiprocessing.Pool() as pool:
pool.map(fill_arr, fill_arr_argDicts)
return shared_array
def write_shp(path_out, shapely_geom, prj=None, attrDict=None):
shapely_geom = [shapely_geom] if not isinstance(shapely_geom, list) else shapely_geom
attrDict = [attrDict] if not isinstance(attrDict, list) else attrDict
# print(len(shapely_geom))
# print(len(attrDict))
assert len(shapely_geom) == len(attrDict), "'shapely_geom' and 'attrDict' must have the same length."
assert os.path.exists(os.path.dirname(path_out)), 'Directory %s does not exist.' % os.path.dirname(path_out)
print('Writing %s ...' % path_out)
if os.path.exists(path_out):
os.remove(path_out)
ds = ogr.GetDriverByName("Esri Shapefile").CreateDataSource(path_out)
if prj is not None:
prj = prj if not isinstance(prj, int) else EPSG2WKT(prj)
srs = osr.SpatialReference()
srs.ImportFromWkt(prj)
else:
srs = None
geom_type = list(set([gm.type for gm in shapely_geom]))
assert len(geom_type) == 1, 'All shapely geometries must belong to the same type. Got %s.' % geom_type
layer = \
ds.CreateLayer('', srs, ogr.wkbPoint) if geom_type[0] == 'Point' else\
ds.CreateLayer('', srs, ogr.wkbLineString) if geom_type[0] == 'LineString' else \
ds.CreateLayer('', srs, ogr.wkbPolygon) if geom_type[0] == 'Polygon' else None # FIXME
if isinstance(attrDict[0], dict):
for attr in attrDict[0].keys():
assert len(attr) <= 10, "ogr does not support fieldnames longer than 10 digits. '%s' is too long" % attr
DTypeStr = get_dtypeStr(attrDict[0][attr])
FieldType = \
ogr.OFTInteger if DTypeStr.startswith('int') else \
ogr.OFTReal if DTypeStr.startswith('float') else \
ogr.OFTString if DTypeStr.startswith('str') else \
ogr.OFTDateTime if DTypeStr.startswith('date') else None
FieldDefn = ogr.FieldDefn(attr, FieldType)
if DTypeStr.startswith('float'):
FieldDefn.SetPrecision(6)
layer.CreateField(FieldDefn) # Add one attribute
for i in range(len(shapely_geom)):
# Create a new feature (attribute and geometry)
feat = ogr.Feature(layer.GetLayerDefn())
feat.SetGeometry(ogr.CreateGeometryFromWkb(shapely_geom[i].wkb)) # Make a geometry, from Shapely object
list_attr2set = attrDict[0].keys() if isinstance(attrDict[0], dict) else []
for attr in list_attr2set:
val = attrDict[i][attr]
DTypeStr = get_dtypeStr(val)
val = int(val) if DTypeStr.startswith('int') else float(val) if DTypeStr.startswith('float') else \
str(val) if DTypeStr.startswith('str') else val
feat.SetField(attr, val)
layer.CreateFeature(feat)
feat.Destroy()
# Save and close everything
del ds, layer
def write_numpy_to_image(array, path_out, outFmt='GTIFF', gt=None, prj=None):
rows, cols, bands = list(array.shape) + [1] if len(array.shape) == 2 else array.shape
gdal_dtype = gdal.GetDataTypeByName(convertGdalNumpyDataType(array.dtype))
outDs = gdal.GetDriverByName(outFmt).Create(path_out, cols, rows, bands, gdal_dtype)
for b in range(bands):
band = outDs.GetRasterBand(b + 1)
arr2write = array if len(array.shape) == 2 else array[:, :, b]
band.WriteArray(arr2write)
del band
if gt:
outDs.SetGeoTransform(gt)
if prj: