Commit 892abd91 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Initial commit

parents
# Created by .ignore support plugin (hsz.mobi)
.idea/
\ No newline at end of file
from . import io
from . import map
from . import numeric
from . import processing
__all__=['io',
'map',
'numeric',
'processing']
__version__ = '20160927.01'
__author__='Daniel Scheffler'
\ No newline at end of file
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
from . import raster
__all__=['raster']
\ No newline at end of file
This diff is collapsed.
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
from . import GeoArray
__all__=['GeoArray']
\ No newline at end of file
from . import raster
from . import vector
from . import coord_calc
from . import coord_trafo
from . import coord_grid
from . import map_info
from . import projection
__all__=['raster',
'vector',
'coord_calc',
'coord_trafo',
'coord_grid',
'map_info',
'projection']
__author__='Daniel Scheffler'
\ No newline at end of file
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
import numpy as np
# custom
from shapely.geometry import Polygon
from .vector import topology as TOPO
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 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(TOPO.find_line_intersection_point(line_N, line_W))),
list(reversed(TOPO.find_line_intersection_point(line_N, line_O))),
list(reversed(TOPO.find_line_intersection_point(line_S, line_W))),
list(reversed(TOPO.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:
# envi.save_image(os.path.abspath('./testing/out/__faulty_mask_1bit.hdr'), self.mask_1bit, dtype= 'uint8', interleave='bsq', ext='.bsq',force=True)
# print(UL,UR,LL,LR)
return out
\ No newline at end of file
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
import numpy as np
from ..numeric.vector import find_nearest
def snap_bounds_to_pixGrid(bounds, gt, roundAlg='auto'):
in_xmin, in_ymin, in_xmax, in_ymax = bounds
xgrid = np.arange(gt[0], gt[0] + 2 * gt[1], gt[1])
ygrid = np.arange(gt[3], gt[3] + 2 * abs(gt[5]), abs(gt[5]))
xmin = find_nearest(xgrid, in_xmin, roundAlg, extrapolate=True)
ymax = find_nearest(ygrid, in_ymax, roundAlg, extrapolate=True)
xmax = find_nearest(xgrid, in_xmax, roundAlg, extrapolate=True)
ymin = find_nearest(ygrid, in_ymin, roundAlg, extrapolate=True)
return xmin, ymin, xmax, ymax
def is_coord_grid_equal(gt,xgrid,ygrid):
ULx,xgsd,holder,ULy,holder,ygsd = gt
if is_point_on_grid((ULx,ULy),xgrid,ygrid) and xgsd==abs(xgrid[1]-xgrid[0]) and abs(ygsd)==abs(ygrid[1]-ygrid[0]):
return True
else:
return False
def is_point_on_grid(pointXY,xgrid,ygrid):
# type: (tuple,np.array,np.array)
"""Checks if a given point is exactly on the given coordinate grid.
:param pointXY: (X,Y) coordinates of the point to check
:param xgrid: numpy array defining the coordinate grid in x-direction
:param ygrid: numpy array defining the coordinate grid in y-direction
"""
if find_nearest(xgrid,pointXY[0],extrapolate=True)==pointXY[0] and \
find_nearest(ygrid,pointXY[1],extrapolate=True)==pointXY[1]:
return True
else:
return False
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
import pyproj
try:
from osgeo import ogr
from osgeo import osr
from osgeo import gdal
from osgeo import gdalnumeric
from osgeo.gdalconst import *
except ImportError:
import ogr
import osr
import gdal
import gdalnumeric
from gdalconst import *
from .projection import get_proj4info
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 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 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 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
\ No newline at end of file
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
try:
from osgeo import osr
except ImportError:
import osr
from .coord_trafo import transform_utm_to_wgs84
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
\ No newline at end of file
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
import re
import os
# custom
try:
from osgeo import ogr
from osgeo import osr
from osgeo import gdal
from osgeo import gdalnumeric
from osgeo.gdalconst import *
except ImportError:
import ogr
import osr
import gdal
import gdalnumeric
from gdalconst import *
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 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
\ No newline at end of file
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"