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

moved a couple of functions from GEOPROCESSING to external package...

moved a couple of functions from GEOPROCESSING to external package 'py_tools_ds' and added direct imports within respective modules
GEOP:
- the following functions moved to 'py_tools_ds':
    - transform_utm_to_wgs84
    - transform_wgs84_to_utm
    - transform_any_prj
    - reproject_shapelyPoly
    - lonlat_to_pixel
    - latLonToPixel
    - pixelToLatLon
    - pixelToMapYX
    - mapXY2imXY
    - imXY2mapXY
    - isProjectedOrGeographic
    - EPSG2Proj4
    - EPSG2WKT
    - WKT2EPSG
    - get_UTMzone
    - geotransform2mapinfo
    - mapinfo2geotransform
    - get_corner_coordinates
    - get_prjLonLat
    - get_proj4info
    - proj4_to_dict
    - prj_equal
    - corner_coord_to_minmax
    - get_footprint_polygon
    - get_overlap_polygon
    - find_line_intersection_point
    - calc_FullDataset_corner_positions
    - is_point_on_grid
    - is_coord_grid_equal
    - snap_bounds_to_pixGrid
IO:
- deleted deprecated module envifilehandling
- renamed envifilehandling_BD to envifilehandling
- deleted deprecated module sysenvironment
parent be17e860
......@@ -27,7 +27,6 @@ import math
import subprocess
import builtins
import time
import warnings
import collections
from scipy.interpolate import RegularGridInterpolator
......@@ -48,14 +47,14 @@ except ImportError:
import pyproj
from pyorbital import astronomy
import ephem
from spectral.io import envi
from shapely.geometry import Polygon
from shapely.geometry import shape
from shapely.geometry import MultiPoint
from numba import jit
from ..io import envifilehandling_BD as ef
from ..io import envifilehandling as ef
from ..misc import helper_functions as HLP_F
from py_tools_ds.ptds.geo.coord_grid import snap_bounds_to_pixGrid
from py_tools_ds.ptds.geo.coord_trafo import transform_utm_to_wgs84, transform_wgs84_to_utm, mapXY2imXY, imXY2mapXY
from py_tools_ds.ptds.geo.projection import get_UTMzone, EPSG2WKT, isProjectedOrGeographic
from py_tools_ds.ptds.geo.raster.reproject import warp_ndarray
job, usecase = builtins.GMS_config.job, builtins.GMS_config.usecase # read from builtins (set by process_controller)
......@@ -2702,59 +2701,6 @@ def Intersect_ext_hdr(hdr_l, ref_ext, v=0):
return hdr_inters
def transform_utm_to_wgs84(easting, northing, zone, south=False):
# ''' returns lon, lat, altitude '''
# utm_coordinate_system = osr.SpatialReference()
# utm_coordinate_system.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon
# is_northern = northingcalc_FullDataset_corner_positions > 0
# utm_coordinate_system.SetUTM(zone, is_northern)
#
# wgs84_coordinate_system = utm_coordinate_system.CloneGeogCS() # Clone ONLY the geographic coordinate system
#
# # create transform component
# utm_to_wgs84_geo_transform = osr.CoordinateTransformation(utm_coordinate_system, wgs84_coordinate_system) # (, )
# print('hier2')
# print(utm_to_wgs84_geo_transform.TransformPoint(easting, northing, 0))
# return utm_to_wgs84_geo_transform.TransformPoint(easting, northing, 0) #[Lon,Lat,meters above sea level]
UTM = pyproj.Proj(proj='utm', zone=abs(zone), ellps='WGS84',south=(zone<0 or south))
return UTM(easting, northing, inverse=True)
def transform_wgs84_to_utm(lon, lat, zone=None):
""" returns easting, northing, altitude. If zone is not set it is derived automatically.
:param lon:
:param lat:
:param zone:
"""
assert -180.<=lon<=180. and -90.<=lat<=90., '''GEOPROCESSING error: Input coordinates for transform_wgs84_to_utm()
out of range. Got %s/%s.''' %(lon,lat)
get_utm_zone = lambda longitude: int(1 + (longitude + 180.0) / 6.0)
# def is_northern(latitude):
# """
# Determines if given latitude is a northern for UTM
# """
# if (latitude < 0.0):
# return 0
# else:
# return 1
#
# utm_coordinate_system = osr.SpatialReference()
# utm_coordinate_system.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon
# utm_coordinate_system.SetUTM(get_utm_zone(lon), is_northern(lat))
#
# wgs84_coordinate_system = utm_coordinate_system.CloneGeogCS() # Clone ONLY the geographic coordinate system
#
# # create transform component
# wgs84_to_utm_geo_transform = osr.CoordinateTransformation(wgs84_coordinate_system, utm_coordinate_system) # (, )
# return wgs84_to_utm_geo_transform.TransformPoint(lon, lat, 0) # [RW,HW,meters above sea level]
UTM = pyproj.Proj(proj='utm', zone=get_utm_zone(lon) if zone is None else zone, ellps='WGS84')
return UTM(lon, lat)
def transform_any_prj(prj_src,prj_tgt,x,y):
prj_src = pyproj.Proj(get_proj4info(proj=prj_src))
prj_tgt = pyproj.Proj(get_proj4info(proj=prj_tgt))
x,y = pyproj.transform(prj_src,prj_tgt,x,y)
return x,y
def get_point_bounds(points):
mp_points = MultiPoint(points)
......@@ -2762,25 +2708,6 @@ def get_point_bounds(points):
return min_lon, min_lat, max_lon, max_lat
def reproject_shapelyPoly(shapelyPoly, tgt_prj):
# type: (shapely.Polygon, str) -> shapely.Polygon
"""Repojects a shapely polygon that has LonLat coordinates into the given target projection.
:param shapelyPoly: <shapely.Polygon> the shapely polygon to be reprojected (must have lonlat coordinates)
:param tgt_prj: <str> WKT projection string (like GDAL projection)
"""
# TODO nicht sicher, ob wirklich nur lonlat funktioniert
get_coordsArr = lambda shpPoly: np.swapaxes(np.array(shpPoly.exterior.coords.xy), 0, 1)
coordsArr = get_coordsArr(shapelyPoly)
coordsArrTgtPrj = np.zeros_like(coordsArr)
proj4dict = proj4_to_dict(get_proj4info(proj=tgt_prj))
ReProj = pyproj.Proj(proj=proj4dict['proj'], zone=proj4dict['zone'], ellps=proj4dict['datum'])
x, y = ReProj(coordsArr[:, 0], coordsArr[:, 1], inverse=False)
coordsArrTgtPrj[:, 0] = x
coordsArrTgtPrj[:, 1] = y
return Polygon(coordsArrTgtPrj)
def zoom_2Darray_to_shapeFullArr(arr2D,shapeFullArr,subset=None,method='linear'):
assert method in ['linear','nearest']
rS,rE,cS,cE = list(get_subsetProps_from_subsetArg(shapeFullArr,subset).values())[3:7]
......@@ -2791,6 +2718,7 @@ def zoom_2Darray_to_shapeFullArr(arr2D,shapeFullArr,subset=None,method='linear')
out_rows_grid, out_cols_grid = np.meshgrid(range(rS,rE),range(cS,cE),indexing='ij')
return rgi(np.dstack([out_rows_grid,out_cols_grid]))
def adjust_acquisArrProv_to_shapeFullArr(arrProv,shapeFullArr,subset=None,bandwise=False):
assert isinstance(arrProv,dict), 'Expected a dictionary with band names (from LayerbandsAssignment) as keys and ' \
'arrays as values. Got %s.' %type(arrProv)
......@@ -2806,347 +2734,13 @@ def adjust_acquisArrProv_to_shapeFullArr(arrProv,shapeFullArr,subset=None,bandwi
interpolated = zoom_2Darray_to_shapeFullArr(arr2interp,shapeFullArr,subset).astype(np.float32)
return interpolated
def get_raster_size(minx, miny, maxx, maxy, cell_width, cell_height):
"""Determine the number of rows/columns given the bounds of the point data and the desired cell size"""
cols = int((maxx - minx) / cell_width)
rows = int((maxy - miny) / abs(cell_height))
return cols, rows
def lonlat_to_pixel(lon, lat, inverse_geo_transform):
"""Translates the given lon, lat to the grid pixel coordinates in data array (zero start)"""
# transform to utm
utm_x, utm_y, utm_z = transform_wgs84_to_utm(lon, lat)
# apply inverse geo tranform
pixel_x, pixel_y = gdal.ApplyGeoTransform(inverse_geo_transform, utm_x, utm_y)
pixel_x = int(pixel_x) - 1 # adjust to 0 start for array
pixel_y = int(pixel_y) - 1 # adjust to 0 start for array
return pixel_x, abs(pixel_y) # y pixel is likly a negative value given geo_transform
def latLonToPixel(latLonPairs, path_im=None, geotransform=None, projection=None):
"""The following method translates given latitude/longitude pairs into pixel locations on a given GEOTIF
This method does not take into account pixel size and assumes a high enough
image resolution for pixel size to be insignificant
:param latLonPairs: The decimal lat/lon pairings to be translated in the form [[lat1,lon1],[lat2,lon2]]
:param path_im: The file location of the input image
:param geotransform: GDAL GeoTransform
:param projection: GDAL Projection
:return: The pixel translation of the lat/lon pairings in the form [[x1,y1],[x2,y2]]
"""
assert path_im or geotransform and projection, \
"GEOP.latLonToPixel: Missing argument! Please provide either 'path_im' or 'geotransform' AND 'projection'."
if geotransform and projection:
gt,proj = geotransform, projection
else:
ds = gdal.Open(path_im)
gt, proj = ds.GetGeoTransform(), ds.GetProjection()
# Load the image dataset
# Create a spatial reference object for the dataset
srs = osr.SpatialReference()
srs.ImportFromWkt(proj)
# Set up the coordinate transformation object
srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srsLatLong, srs)
# Go through all the point pairs and translate them to latitude/longitude pairings
pixelPairs = []
for point in latLonPairs:
# Change the point locations into the GeoTransform space
(point[1], point[0], holder) = ct.TransformPoint(point[1], point[0])
# Translate the x and y coordinates into pixel values
x = (point[1] - gt[0]) / gt[1]
y = (point[0] - gt[3]) / gt[5]
# Add the point to our return array
pixelPairs.append([int(x), int(y)])
return pixelPairs
def pixelToLatLon(pixelPairs, path_im=None, geotransform=None, projection=None):
"""The following method translates given pixel locations into latitude/longitude locations on a given GEOTIF
This method does not take into account pixel size and assumes a high enough
image resolution for pixel size to be insignificant
:param pixelPairs: Image coordinate pairs to be translated in the form [[x1,y1],[x2,y2]]
:param path_im: The file location of the input image
:param geotransform: GDAL GeoTransform
:param projection: GDAL Projection
:returns: The lat/lon translation of the pixel pairings in the form [[lat1,lon1],[lat2,lon2]]
"""
assert path_im or geotransform and projection, \
"GEOP.pixelToLatLon: Missing argument! Please provide either 'path_im' or 'geotransform' AND 'projection'."
if geotransform and projection:
gt,proj = geotransform, projection
else:
ds = gdal.Open(path_im)
gt, proj = ds.GetGeoTransform(), ds.GetProjection()
# Create a spatial reference object for the dataset
srs = osr.SpatialReference()
srs.ImportFromWkt(proj)
# Set up the coordinate transformation object
srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srs, srsLatLong)
# Go through all the point pairs and translate them to pixel pairings
latLonPairs = []
for point in pixelPairs:
# Translate the pixel pairs into untranslated points
ulon = point[0] * gt[1] + gt[0]
ulat = point[1] * gt[5] + gt[3]
# Transform the points to the space
(lon, lat, holder) = ct.TransformPoint(ulon, ulat)
# Add the point to our return array
latLonPairs.append([lat, lon])
return latLonPairs
def pixelToMapYX(pixelCoords, path_im=None, geotransform=None, projection=None):
"""Translates given pixel locations into map coordinates of the given image.
:param pixelCoords: The pixel coordinates to be translated in the form [x,y].
A list of coordinate pairs is possible: [[x1,y1],[x2,y2]].
:param path_im: absolute path of the image to be used to get geotransform and projection information
:param geotransform: GDAL geotransform
:param projection: GDAL projection
:return: The lat/lon or UTM translation result of the pixel coordinates
in the form [lat,lon] or [[lat1,lon1],[lat2,lon2]].
"""
assert path_im or geotransform and projection, \
"GEOP.pixelToMapYX: Missing argument! Please provide either 'path_im' or 'geotransform' AND 'projection'."
if geotransform and projection:
gt,proj = geotransform, projection
else:
ds = gdal.Open(path_im)
gt, proj = ds.GetGeoTransform(), ds.GetProjection()
srs = osr.SpatialReference()
srs.ImportFromWkt(proj)
# Set up the coordinate transformation object
ct = osr.CoordinateTransformation(srs, srs)
mapYmapXPairs = []
pixelCoords = [pixelCoords] if not isinstance(pixelCoords[0],list) else pixelCoords
for point in pixelCoords:
# Translate the pixel pairs into untranslated points
u_mapX = point[0] * gt[1] + gt[0]
u_mapY = point[1] * gt[5] + gt[3]
# Transform the points to the space
(mapX, mapY, holder) = ct.TransformPoint(u_mapX, u_mapY)
# Add the point to our return array
mapYmapXPairs.append([mapY, mapX])
return mapYmapXPairs
def mapXY2imXY(mapXY, gt):
# type: (tuple,list) -> tuple
"""Translates given map coordinates into pixel locations according to the given image geotransform.
:param mapXY: <tuple> The map coordinates to be translated in the form (x,y).
:param gt: <list> GDAL geotransform
:returns: <tuple> image coordinate tuple X/Y (column, row)
"""
return (mapXY[0]-gt[0])/gt[1],(mapXY[1]-gt[3])/gt[5]
def imXY2mapXY(imXY, gt):
# type: (tuple,list) -> tuple
"""Translates given pixel locations into map coordinates according to the given image geotransform.
:param imXY: <tuple> The image coordinates to be translated in the form (x,y).
:param gt: <list> GDAL geotransform
:returns: <tuple> map coordinate tuple X/Y (mapX, mapY)
"""
return (gt[0] + imXY[0] * abs(gt[1])), (gt[3] - imXY[1] * abs(gt[5]))
def isProjectedOrGeographic(prj):
""":param prj: accepts EPSG, Proj4 and WKT projections"""
srs = osr.SpatialReference()
if prj.startswith('EPSG:'):
srs.ImportFromEPSG(int(prj.split(':')[1]))
elif prj.startswith('+proj='):
srs.ImportFromProj4(prj)
elif prj.startswith('GEOGCS') or prj.startswith('PROJCS'):
srs.ImportFromWkt(prj)
else:
raise Exception('Unknown input projection.')
return 'projected' if srs.IsProjected() else 'geographic' if srs.IsGeographic() else None
def EPSG2Proj4(EPSG_code):
srs = osr.SpatialReference()
srs.ImportFromEPSG(EPSG_code)
return srs.ExportToProj4()
def EPSG2WKT(EPSG_code):
srs = osr.SpatialReference()
srs.ImportFromEPSG(EPSG_code)
return srs.ExportToWkt()
def WKT2EPSG(wkt, epsg=os.environ['GDAL_DATA'].replace('/gdal','/proj/epsg')):
""" Transform a WKT string to an EPSG code
:param wkt: WKT definition, e.g. the result of gdal_ds.GetProjection()
:param epsg: the proj.4 epsg file (defaults to '/usr/local/share/proj/epsg')
:returns: EPSG code
http://gis.stackexchange.com/questions/20298/is-it-possible-to-get-the-epsg-value-from-an-osr-spatialreference-class-using-th
"""
code = None #default
p_in = osr.SpatialReference()
s = p_in.ImportFromWkt(wkt)
if s == 5: # invalid WKT
raise Exception('Invalid WKT.')
if p_in.IsLocal():
raise Exception('The given WKT is a local coordinate system.')
cstype = 'GEOGCS' if p_in.IsGeographic() else 'PROJCS'
an = p_in.GetAuthorityName(cstype)
assert an in [None,'EPSG'], "No EPSG code found. Found %s instead." %an
ac = p_in.GetAuthorityCode(cstype)
if ac is None: # try brute force approach by grokking proj epsg definition file
p_out = p_in.ExportToProj4()
if p_out:
with open(epsg) as f:
for line in f:
if line.find(p_out) != -1:
m = re.search('<(\\d+)>', line)
if m:
code = m.group(1)
break
code = code if code else None # match or no match
else:
code = ac
return int(code)
def get_UTMzone(ds=None,prj=None):
assert ds or prj, "Specify at least one of the arguments 'ds' or 'prj'"
if isProjectedOrGeographic(prj)=='projected':
srs = osr.SpatialReference()
srs.ImportFromWkt(ds.GetProjection() if ds else prj)
return srs.GetUTMZone()
else:
return None
def geotransform2mapinfo(gt,prj):
"""Builds an ENVI map info from given GDAL GeoTransform and Projection (compatible with UTM and LonLat projections).
:param gt: GDAL GeoTransform, e.g. (249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0)
:param prj: GDAL Projection - WKT Format
:returns: ENVI map info, e.g. [ UTM , 1 , 1 , 256785.0 , 4572015.0 , 30.0 , 30.0 , 43 , North , WGS-84 ]
:rtype: list
"""
if gt[2]!=0 or gt[4]!=0: # TODO
raise NotImplementedError('Currently rotated datasets are not supported.')
srs = osr.SpatialReference()
srs.ImportFromWkt(prj)
Proj4 = [i[1:] for i in srs.ExportToProj4().split()]
Proj4_proj = [v.split('=')[1] for i,v in enumerate(Proj4) if '=' in v and v.split('=')[0]=='proj'][0]
Proj4_ellps = [v.split('=')[1] for i,v in enumerate(Proj4) if '=' in v and v.split('=')[0] in ['ellps','datum']][0]
proj = 'Geographic Lat/Lon' if Proj4_proj=='longlat' else 'UTM' if Proj4_proj=='utm' else Proj4_proj
ellps = 'WGS-84' if Proj4_ellps=='WGS84' else Proj4_ellps
UL_X, UL_Y, GSD_X, GSD_Y = gt[0], gt[3], gt[1], gt[5]
utm2wgs84 = lambda utmX,utmY: transform_utm_to_wgs84(utmX,utmY,srs.GetUTMZone())
is_UTM_North_South = lambda LonLat: 'North' if LonLat[1] >= 0. else 'South'
map_info = [proj,1,1,UL_X,UL_Y,GSD_X,abs(GSD_Y),ellps] if proj!='UTM' else \
['UTM',1,1,UL_X,UL_Y,GSD_X,abs(GSD_Y),srs.GetUTMZone(),is_UTM_North_South(utm2wgs84(UL_X,UL_Y)),ellps]
return map_info
def mapinfo2geotransform(map_info):
# type: (list) -> list
"""Builds GDAL GeoTransform tuple from an ENVI map info.
:param map_info: ENVI map info (list), e.g. ['UTM', 1, 1, 192585.0, 5379315.0, 30.0, 30.0, 41, 'North', 'WGS-84']
:returns: GDAL GeoTransform, e.g. [249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0]
:rtype: list
"""
# 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:
ULmapX, ULmapY = float(map_info[3]),float(map_info[4])
else:
ULmapX = float(map_info[3])-(float(map_info[1])*float(map_info[5])-float(map_info[5]))
ULmapY = float(map_info[4])+(float(map_info[2])*float(map_info[6])-float(map_info[6]))
assert isinstance(map_info,list) and len(map_info)==10,\
"The map info argument has to be a list of length 10. Got %s." %map_info
return [ULmapX,float(map_info[5]),0.,ULmapY,0.,-float(map_info[6])]
def get_corner_coordinates(gdal_ds=None, gt=None, cols=None, rows=None):
"""Returns (ULxy, LLxy, LRxy, URxy) in the same coordinate units like the given geotransform."""
assert gdal_ds or (gt and cols and rows), \
"GEOP.get_corner_coordinates: Missing argument! Please provide either 'gdal_ds' or 'gt', 'cols' AND 'rows'."
gdal_ds_GT = gdal_ds.GetGeoTransform() if gdal_ds else gt
ext=[]
xarr=[0,gdal_ds.RasterXSize if gdal_ds else cols]
yarr=[0,gdal_ds.RasterYSize if gdal_ds else rows]
for px in xarr:
for py in yarr:
x=gdal_ds_GT[0]+(px*gdal_ds_GT[1])+(py*gdal_ds_GT[2])
y=gdal_ds_GT[3]+(px*gdal_ds_GT[4])+(py*gdal_ds_GT[5])
ext.append([x,y])
yarr.reverse()
gdal_ds_GT = None
return ext
def get_prjLonLat(fmt='wkt'):
# type: (str) -> Any
"""Returns standard geographic projection (EPSG 4326) in the WKT or PROJ4 format.
:param fmt: <str> target format - 'WKT' or 'PROJ4'
"""
assert re.search('wkt',fmt,re.I) or re.search('Proj4',fmt,re.I), 'unsupported output format'
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
return srs.ExportToWkt() if re.search('wkt',fmt,re.I) else srs.ExportToProj4()
def get_proj4info(ds=None,proj=None):
# type: (gdal.Dataset,str) -> str
"""Returns PROJ4 formatted projection info for the given gdal.Dataset or projection respectivly,
e.g. '+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs '
:param ds: <gdal.Dataset> the gdal dataset to get PROJ4 info for
:param proj: <str> the projection to get PROJ4 formatted info for
"""
assert ds or proj, "Specify at least one of the arguments 'ds' or 'proj'"
srs = osr.SpatialReference()
srs.ImportFromWkt(ds.GetProjection() if ds else proj)
return srs.ExportToProj4()
def proj4_to_dict(proj4):
#type: (str) -> dict
"""Converts a PROJ4-like string into a dictionary.
:param proj4: <str> the PROJ4-like string
"""
items = [item for item in proj4.replace('+','').split(' ') if '=' in item]
return {k:v for k,v in [kv.split('=') for kv in items]}
def prj_equal(prj1, prj2):
#type: (str,str) -> bool
"""Checks if the given two projections are equal."""
return get_proj4info(proj=prj1)==get_proj4info(proj=prj2)
def corner_coord_to_minmax(corner_coords):
# type: (list) -> (float,float,float,float)
"""Converts [[x1,y1],[x2,y2],[]...] to (xmin,xmax,ymin,ymax)
:param corner_coords: list of coordinates like [[x1,y1],[x2,y2],[]...]
"""
x_vals = [i[0] for i in corner_coords]
y_vals = [i[1] for i in corner_coords]
xmin,xmax,ymin,ymax = min(x_vals),max(x_vals),min(y_vals),max(y_vals)
return xmin,xmax,ymin,ymax
def get_footprint_polygon(CornerLonLat):
""" 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
"""
outpoly = Polygon(CornerLonLat)
assert outpoly. is_valid, 'The given coordinates result in an invalid polygon. Check coordinate order.'
return outpoly
def get_overlap_polygon(poly1, poly2):
""" 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 poly2: second shapely.Polygon() object
:return: overlap polygon as shapely.Polygon() object
:return: overlap percentage as float value [%]
:return: area of overlap polygon
"""
overlap_poly = poly1.intersection(poly2)
if not overlap_poly.is_empty:
overlap_percentage = 100 * shape(overlap_poly).area / shape(poly2).area
return {'overlap poly':overlap_poly, 'overlap percentage':overlap_percentage, 'overlap area':overlap_poly.area}
else:
return {'overlap poly':None, 'overlap percentage':0, 'overlap area':0}
def DN2Rad(ndarray, offsets, gains, inFill=None, inZero=None, inSaturated=None, cutNeg=True):
#type: (np.ndarray,list,list,int,int,int,bool) -> np.ndarray
......@@ -3324,120 +2918,6 @@ def DN2DegreesCelsius_fastforward(ndarray,offsets, gains, K1, K2, emissivity=0.9
return degCel
def find_line_intersection_point(line1, line2):
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])
def det(a, b):
return a[0] * b[1] - a[1] * b[0]
div = det(xdiff, ydiff)
if div == 0:
raise Exception('lines do not intersect')
d = (det(*line1), det(*line2))
x = det(d, xdiff) / div
y = det(d, ydiff) / div
return x, y
def calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=True, algorithm='shapely'):
"""
Calculates the image coordinates of the true data corners for rotated datasets.
ONLY usable for entire images - no tiles!
NOTE: Algorithm calculates the corner coordinates of the convex hull of the given mask. Since the convex hull not
always reflects all of the true corner coordinates the result has a limitation in this regard.
:param mask_1bit: 2D-numpy array 1bit
:param assert_four_corners:
:param algorithm: <str> 'shapely' or 'numpy'
:return: [UL, UR, LL, LR] as [(ULrow,ULcol),(URrow,URcol),...]
"""
rows, cols = mask_1bit.shape[:2]
if sum([mask_1bit[0,0], mask_1bit[0,cols-1], mask_1bit[rows-1,0], mask_1bit[rows-1,cols-1]]) == 4:
# if mask contains pixel value 1 in all of the four corners: return outer corner coords
out = [(0,0),(0,cols-1),(rows-1,0),(rows-1,cols-1)] # UL, UR, LL, LR
else:
assert algorithm in ['shapely','numpy'], "Algorithm '%' does not exist. Choose between 'shapely' and 'numpy'."
if algorithm == 'shapely':
# get outline
topbottom = np.empty((1, 2 * mask_1bit.shape[1]), dtype=np.uint16)
topbottom[0, 0:mask_1bit.shape[1]] = np.argmax(mask_1bit, axis=0)
topbottom[0, mask_1bit.shape[1]:] = (mask_1bit.shape[0] - 1) - np.argmax(np.flipud(mask_1bit), axis=0)
mask = np.tile(np.any(mask_1bit, axis=0), (2,))
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]
if assert_four_corners or len(unique_corn_YX)==4:
# 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)
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()]
out = [getClosest(tgtYX,unique_corn_YX) for tgtYX in [(0,0),(0,cols-1),(rows-1,0),(rows-1,cols-1)]]
else:
out = unique_corn_YX
else: # alg='numpy'
cols_containing_data = np.any(mask_1bit, axis=0)
rows_containing_data = np.any(mask_1bit, axis=1)
first_dataCol = list(cols_containing_data).index(True)
last_dataCol = cols - (list(reversed(cols_containing_data)).index(True) + 1)
first_dataRow = list(rows_containing_data).index(True)
last_dataRow = rows - (list(reversed(rows_containing_data)).index(True) + 1)
StartStopRows_in_first_dataCol = [list(mask_1bit[:, first_dataCol]).index(True),
(rows - list(reversed(mask_1bit[:, first_dataCol])).index(True) + 1)]
StartStopRows_in_last_dataCol = [list(mask_1bit[:, last_dataCol]).index(True),
(rows - list(reversed(mask_1bit[:, last_dataCol])).index(True) + 1)]
StartStopCols_in_first_dataRow = [list(mask_1bit[first_dataRow, :]).index(True),
(cols - list(reversed(mask_1bit[first_dataRow, :])).index(True) + 1)]
StartStopCols_in_last_dataRow = [list(mask_1bit[last_dataRow, :]).index(True),
(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]]:
''' 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]))
line_S = ((first_dataCol, StartStopRows_in_first_dataCol[1]),
(StartStopCols_in_last_dataRow[0], last_dataRow))
line_W = ((StartStopCols_in_first_dataRow[0], first_dataRow),
(first_dataCol, StartStopRows_in_first_dataCol[0]))
line_O = ((last_dataCol, StartStopRows_in_last_dataCol[1]),
(StartStopCols_in_last_dataRow[1], last_dataRow))
corners = [list(reversed(find_line_intersection_point(line_N, line_W))),
list(reversed(find_line_intersection_point(line_N, line_O))),
list(reversed(find_line_intersection_point(line_S, line_W))),
list(reversed(find_line_intersection_point(line_S, line_O)))]
else:
dataRow_in_first_dataCol = np.mean(StartStopRows_in_first_dataCol)
dataRow_in_last_dataCol = np.mean(StartStopRows_in_last_dataCol)
dataCol_in_first_dataRow = np.mean(StartStopCols_in_first_dataRow)
dataCol_in_last_dataRow = np.mean(StartStopCols_in_last_dataRow)
corners = [(first_dataRow, dataCol_in_first_dataRow), (dataRow_in_last_dataCol, last_dataCol),
(last_dataRow, dataCol_in_last_dataRow), (dataRow_in_first_dataCol, first_dataCol)]
getClosest = lambda refR, refC: \
corners[np.array([np.sqrt((refR - c[0]) ** 2 + (refC - c[1]) ** 2) for c in corners]).argmin()]
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: