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

added functions to visualize GeoArray as map; improved GeoArray.show(); some new functions

geo.raster.reproject:
- warp_ndarray():
    - added keywords polynomialOrder,  transformerOptions, warpOptions

geo.vector.geometry:
- boxObj:
    - added get_coordArray_MapXY(): a function to quickly get a geolocation array

geo.coord_trafo:
- transform_any_prj(): added docstring
- added transform_coordArray()
- added reproject_shapelyGeometry()
- reproject_shapelyPoly() is now deprecated

geo.projection:
- get_proj4info() now also accepts EPSG codes

io.raster.GeoArray:
- GeoArray:
    - revised default values of geotransform and projection properties
    - added some validation to geotransform.setter
    - added property epsg
    - revised __getattr__
    - added _get_plottable_image()
    - revised show(): much faster, better colormap
    - added show_map(): a function for quickly showing a map of the GeoArray
    - added show_map_utm() - not yet working
parent 12889c79
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler" __author__ = "Daniel Scheffler"
from functools import partial
import warnings
import pyproj import pyproj
import numpy as np import numpy as np
from shapely.geometry import Polygon from shapely.geometry import Polygon
from shapely.ops import transform
try: try:
from osgeo import ogr from osgeo import ogr
...@@ -72,13 +76,56 @@ def transform_wgs84_to_utm(lon, lat, zone=None): ...@@ -72,13 +76,56 @@ def transform_wgs84_to_utm(lon, lat, zone=None):
return UTM(lon, lat) return UTM(lon, lat)
def transform_any_prj(prj_src,prj_tgt,x,y): def c(prj_src,prj_tgt,x,y):
"""Transforms X/Y data from any source projection to any target projection.
:param prj_src: GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)
:param prj_tgt: GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)
:param x:
:param y:
:return:
"""
prj_src = pyproj.Proj(get_proj4info(proj=prj_src)) prj_src = pyproj.Proj(get_proj4info(proj=prj_src))
prj_tgt = pyproj.Proj(get_proj4info(proj=prj_tgt)) prj_tgt = pyproj.Proj(get_proj4info(proj=prj_tgt))
x,y = pyproj.transform(prj_src,prj_tgt,x,y) x,y = pyproj.transform(prj_src,prj_tgt,x,y)
return x,y return x,y
def transform_coordArray(prj_src, prj_tgt, Xarr, Yarr, Zarr=None):
"""Transforms geolocation arrays from one projection into another.
HINT: This function is faster than transform_any_prj but works only for geolocation arrays.
:param prj_src: WKT string
:param prj_tgt: WKT string
:param Xarr: np.ndarray of shape (rows,cols)
:param Yarr: np.ndarray of shape (rows,cols)
:param Zarr: np.ndarray of shape (rows,cols)
:return:
"""
drv = gdal.GetDriverByName('MEM')
geoloc_ds = drv.Create('geoloc',Xarr.shape[1],Xarr.shape[0],3,gdal.GDT_Float64)
geoloc_ds.GetRasterBand(1).WriteArray(Xarr)
geoloc_ds.GetRasterBand(2).WriteArray(Yarr)
if Zarr is not None:
geoloc_ds.GetRasterBand(3).WriteArray(Zarr)
transformer = gdal.Transformer(None, None, ['SRC_SRS=' + prj_src, 'DST_SRS=' + prj_tgt])
status = transformer.TransformGeolocations(geoloc_ds.GetRasterBand(1),geoloc_ds.GetRasterBand(2),
geoloc_ds.GetRasterBand(3))
if status:
raise Exception('Error transforming coordinate array: ' + gdal.GetLastErrorMsg())
Xarr = geoloc_ds.GetRasterBand(1).ReadAsArray()
Yarr = geoloc_ds.GetRasterBand(2).ReadAsArray()
if Zarr is not None:
Zarr = geoloc_ds.GetRasterBand(3).ReadAsArray()
return Xarr, Yarr, Zarr
else:
return Xarr, Yarr
def mapXY2imXY(mapXY, gt): def mapXY2imXY(mapXY, gt):
# type: (tuple,list) -> tuple # type: (tuple,list) -> tuple
"""Translates given geo coordinates into pixel locations according to the given image geotransform. """Translates given geo coordinates into pixel locations according to the given image geotransform.
...@@ -233,12 +280,27 @@ def transform_GCPlist(gcpList, prj_src, prj_tgt): ...@@ -233,12 +280,27 @@ def transform_GCPlist(gcpList, prj_src, prj_tgt):
for GCP in gcpList] for GCP in gcpList]
def reproject_shapelyGeometry(shapelyGeometry, prj_src, prj_tgt):
"""Reprojects any shapely geometry from one projection to another.
:param shapelyGeometry: any shapely geometry instance
:param prj_src: GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)
:param prj_tgt: GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)
"""
project = partial(pyproj.transform, pyproj.Proj(get_proj4info(proj=prj_src)),
pyproj.Proj(get_proj4info(proj=prj_tgt)) )
return transform(project, shapelyGeometry) # apply projection
def reproject_shapelyPoly(shapelyPoly, tgt_prj): def reproject_shapelyPoly(shapelyPoly, tgt_prj):
# type: (shapely.Polygon, str) -> shapely.Polygon # type: (shapely.Polygon, str) -> shapely.Polygon
"""Repojects a shapely polygon that has LonLat coordinates into the given target projection. """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 shapelyPoly: <shapely.Polygon> the shapely polygon to be reprojected (must have lonlat coordinates)
:param tgt_prj: <str> WKT projection string (like GDAL projection) :param tgt_prj: <str> WKT projection string (like GDAL projection)
""" """
warnings.warn('reproject_shapelyPoly() is deprecated. Use reproject_shapelyGeometry() instead.', DeprecationWarning)
# TODO nicht sicher, ob wirklich nur lonlat funktioniert # TODO nicht sicher, ob wirklich nur lonlat funktioniert
get_coordsArr = lambda shpPoly: np.swapaxes(np.array(shpPoly.exterior.coords.xy), 0, 1) get_coordsArr = lambda shpPoly: np.swapaxes(np.array(shpPoly.exterior.coords.xy), 0, 1)
coordsArr = get_coordsArr(shapelyPoly) coordsArr = get_coordsArr(shapelyPoly)
......
...@@ -26,11 +26,18 @@ def get_proj4info(ds=None,proj=None): ...@@ -26,11 +26,18 @@ def get_proj4info(ds=None,proj=None):
"""Returns PROJ4 formatted projection info for the given gdal.Dataset or projection respectivly, """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 ' 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 ds: <gdal.Dataset> the gdal dataset to get PROJ4 info for
:param proj: <str> the projection to get PROJ4 formatted info for :param proj: <str,int> the projection to get PROJ4 formatted info for (WKT or 'epsg:1234' or <EPSG_int>)
""" """
assert ds or proj, "Specify at least one of the arguments 'ds' or 'proj'" assert ds or proj, "Specify at least one of the arguments 'ds' or 'proj'"
srs = osr.SpatialReference() srs = osr.SpatialReference()
srs.ImportFromWkt(ds.GetProjection() if ds else proj) proj = ds.GetProjection() if ds else proj
if isinstance(proj, str) and proj.startswith('epsg:'):
proj = EPSG2WKT(int(proj.split(':')[1]))
elif isinstance(proj,int):
proj = EPSG2WKT(proj)
srs.ImportFromWkt(proj)
return srs.ExportToProj4() return srs.ExportToProj4()
......
...@@ -223,7 +223,8 @@ def warp_GeoArray(geoArr, **kwargs): ...@@ -223,7 +223,8 @@ def warp_GeoArray(geoArr, **kwargs):
def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=(None, None), 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), out_bounds=None, out_bounds_prj=None, out_XYdims = (None,None),
rspAlg='near', in_nodata=None, out_nodata=None, in_alpha=False, rspAlg='near', in_nodata=None, out_nodata=None, in_alpha=False,
out_alpha=False, targetAlignedPixels=False, gcpList=None, options=None, CPUs=1, progress=True, q=False): out_alpha=False, targetAlignedPixels=False, gcpList=None, polynomialOrder=None, options=None,
transformerOptions=None, warpOptions=None, CPUs=1, progress=True, q=False):
""" """
:param ndarray: :param ndarray:
...@@ -248,7 +249,10 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=( ...@@ -248,7 +249,10 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=(
includes the minimum extent. includes the minimum extent.
:param gcpList: <list> list of ground control points in the output coordinate system :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),...] to be used for warping, e.g. [gdal.GCP(mapX,mapY,mapZ,column,row),...]
:param polynomialOrder: <int> order of polynomial GCP interpolation
:param options: <str> additional GDAL options as string, e.g. '-nosrcalpha' or '-order' :param options: <str> additional GDAL options as string, e.g. '-nosrcalpha' or '-order'
:param transformerOptions: <list> list of transformer options, e.g. ['SRC_SRS=invalid']
:param warpOptions: <list> list of warp options, e.g. ['CUTLINE_ALL_TOUCHED=TRUE']
:param CPUs: <int> number of CPUs to use (default: None, which means 'all CPUs available') :param CPUs: <int> number of CPUs to use (default: None, which means 'all CPUs available')
:param progress: <bool> show progress bar (default: True) :param progress: <bool> show progress bar (default: True)
:param q: <bool> quiet mode (default: False) :param q: <bool> quiet mode (default: False)
...@@ -277,9 +281,7 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=( ...@@ -277,9 +281,7 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=(
cropToCutline = False cropToCutline = False
cutlineSQL = 'SELECT * FROM cutline' cutlineSQL = 'SELECT * FROM cutline'
cutlineWhere = '1 = 1' cutlineWhere = '1 = 1'
warpOptions = [] # ['CUTLINE_ALL_TOUCHED=TRUE'] # this is how to implement extra options
callback_data = [0] callback_data = [0]
transformerOptions = ['SRC_SRS=invalid']
rpc = [ rpc = [
"HEIGHT_OFF=1466.05894327379", "HEIGHT_OFF=1466.05894327379",
"HEIGHT_SCALE=144.837606185489", "HEIGHT_SCALE=144.837606185489",
...@@ -298,6 +300,47 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=( ...@@ -298,6 +300,47 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=(
] ]
WarpMemoryLimit = 0 # not sure if this is the correct keyword?? WarpMemoryLimit = 0 # not sure if this is the correct keyword??
""" Create a WarpOptions() object that can be passed to gdal.Warp()
Keyword arguments are :
options --- can be be an array of strings, a string or let empty and filled from other keywords.
format --- output format ("GTiff", etc...)
outputBounds --- output bounds as (minX, minY, maxX, maxY) in target SRS
outputBoundsSRS --- SRS in which output bounds are expressed, in the case they are not expressed in dstSRS
xRes, yRes --- output resolution in target SRS
targetAlignedPixels --- whether to force output bounds to be multiple of output resolution
width --- width of the output raster in pixel
height --- height of the output raster in pixel
srcSRS --- source SRS
dstSRS --- output SRS
srcAlpha --- whether to force the last band of the input dataset to be considered as an alpha band
dstAlpha --- whether to force the creation of an output alpha band
outputType --- output type (gdal.GDT_Byte, etc...)
workingType --- working type (gdal.GDT_Byte, etc...)
warpOptions --- list of warping options
errorThreshold --- error threshold for approximation transformer (in pixels)
warpMemoryLimit --- size of working buffer in bytes
resampleAlg --- resampling mode
creationOptions --- list of creation options
srcNodata --- source nodata value(s)
dstNodata --- output nodata value(s)
multithread --- whether to multithread computation and I/O operations
tps --- whether to use Thin Plate Spline GCP transformer
rpc --- whether to use RPC transformer
geoloc --- whether to use GeoLocation array transformer
polynomialOrder --- order of polynomial GCP interpolation
transformerOptions --- list of transformer options
cutlineDSName --- cutline dataset name
cutlineLayer --- cutline layer name
cutlineWhere --- cutline WHERE clause
cutlineSQL --- cutline SQL statement
cutlineBlend --- cutline blend distance in pixels
cropToCutline --- whether to use cutline extent for output bounds
copyMetadata --- whether to copy source metadata
metadataConflictValue --- metadata data conflict value
setColorInterpretation --- whether to force color interpretation of input bands to output bands
callback --- callback method
callback_data --- user data for callback
"""
# get input dataset (in-MEM) # get input dataset (in-MEM)
...@@ -339,10 +382,11 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=( ...@@ -339,10 +382,11 @@ def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=(
srcAlpha = in_alpha, srcAlpha = in_alpha,
dstAlpha = out_alpha, dstAlpha = out_alpha,
options = options if options else [], options = options if options else [],
#warpOptions = [], warpOptions = warpOptions,
#transformerOptions = [], transformerOptions = transformerOptions,
targetAlignedPixels = targetAlignedPixels, targetAlignedPixels = targetAlignedPixels,
tps = True if gcpList else False, tps = True if gcpList else False,
polynomialOrder = polynomialOrder,
callback = progressBarWarp, callback = progressBarWarp,
errorThreshold = 0.125, # this is needed to get exactly the same output like the console version of GDAL warp errorThreshold = 0.125, # this is needed to get exactly the same output like the console version of GDAL warp
) )
......
...@@ -3,12 +3,14 @@ __author__ = "Daniel Scheffler" ...@@ -3,12 +3,14 @@ __author__ = "Daniel Scheffler"
from shapely.geometry import Polygon, box from shapely.geometry import Polygon, box
import numpy as np
from .conversion import shapelyImPoly_to_shapelyMapPoly, get_boxImXY_from_shapelyPoly, \ from .conversion import shapelyImPoly_to_shapelyMapPoly, get_boxImXY_from_shapelyPoly, \
shapelyBox2BoxYX, round_shapelyPoly_coords shapelyBox2BoxYX, round_shapelyPoly_coords
from ..coord_calc import corner_coord_to_minmax from ..coord_calc import corner_coord_to_minmax
from ..coord_trafo import imYX2mapYX from ..coord_trafo import imYX2mapYX, transform_coordArray
from ..projection import prj_equal
class boxObj(object): class boxObj(object):
...@@ -144,6 +146,22 @@ class boxObj(object): ...@@ -144,6 +146,22 @@ class boxObj(object):
y_is_larger = ymin<b2t_ymin or ymax>b2t_ymax y_is_larger = ymin<b2t_ymin or ymax>b2t_ymax
return x_is_larger,y_is_larger return x_is_larger,y_is_larger
def get_coordArray_MapXY(self, prj=None):
"""Returns two coordinate arrays for X and Y coordinates in the given projection. If no projection is given,
<boxObj>.prj is used.
:param prj: GDAL projection as WKT string
:return:
"""
assert self.gt and self.prj
xmin, xmax, ymin, ymax = self.boundsMap
Xarr, Yarr = np.meshgrid(np.arange(xmin, xmax, abs(self.gt[1])),
np.arange(ymin, ymax, abs(self.gt[5])) )
if prj and not prj_equal(self.prj, prj):
Xarr, Yarr = transform_coordArray(self.prj, prj, Xarr, Yarr)
return Xarr, Yarr
def get_winPoly(wp_imYX, ws, gt, match_grid=0): def get_winPoly(wp_imYX, ws, gt, match_grid=0):
# type: (tuple, tuple, list, bool) -> shapely.Polygon, tuple, tuple # type: (tuple, tuple, list, bool) -> shapely.Polygon, tuple, tuple
......
...@@ -6,6 +6,7 @@ import numpy as np ...@@ -6,6 +6,7 @@ import numpy as np
import os import os
import warnings import warnings
from matplotlib import pyplot as plt from matplotlib import pyplot as plt
from mpl_toolkits.basemap import Basemap
# custom # custom
from shapely.geometry import Polygon, box from shapely.geometry import Polygon, box
...@@ -23,7 +24,7 @@ except ImportError: ...@@ -23,7 +24,7 @@ except ImportError:
from ...geo.coord_calc import get_corner_coordinates, calc_FullDataset_corner_positions from ...geo.coord_calc import get_corner_coordinates, calc_FullDataset_corner_positions
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, WKT2EPSG, EPSG2WKT
from ...geo.vector.topology import get_overlap_polygon, get_footprint_polygon from ...geo.vector.topology import get_overlap_polygon, get_footprint_polygon
from ...geo.vector.geometry import boxObj from ...geo.vector.geometry import boxObj
from ...io.raster.gdal import get_GDAL_ds_inmem from ...io.raster.gdal import get_GDAL_ds_inmem
...@@ -142,11 +143,14 @@ class GeoArray(object): ...@@ -142,11 +143,14 @@ class GeoArray(object):
self.set_gdalDataset_meta() self.set_gdalDataset_meta()
return self._geotransform return self._geotransform
else: else:
raise AttributeError("Attribute 'geotransform' has not been set yet.") return [0,1,0,0,0,-1]
@geotransform.setter @geotransform.setter
def geotransform(self, gt): def geotransform(self, gt):
assert isinstance(gt,(list,tuple)) and len(gt)==6, 'geotransform must be a list with 6 numbers. Got %s.' %gt
for i in gt: assert isinstance(i,(int,float)), "geotransform must contain only numbers. Got '%s'." %i
if self.filePath: if self.filePath:
assert self.geotransform == gt, "Cannot set %s.geotransform to the given value because it does not " \ assert self.geotransform == gt, "Cannot set %s.geotransform to the given value because it does not " \
"match the geotransform from the file on disk." %self.__class__.__name__ "match the geotransform from the file on disk." %self.__class__.__name__
...@@ -172,7 +176,7 @@ class GeoArray(object): ...@@ -172,7 +176,7 @@ class GeoArray(object):
self.set_gdalDataset_meta() self.set_gdalDataset_meta()
return self._projection return self._projection
else: else:
raise AttributeError("Attribute 'projection' has not been set yet.") return ''
@projection.setter @projection.setter
...@@ -184,6 +188,16 @@ class GeoArray(object): ...@@ -184,6 +188,16 @@ class GeoArray(object):
self._projection = prj self._projection = prj
@property
def epsg(self):
return WKT2EPSG(self.projection)
@epsg.setter
def epsg(self, epsg_code):
self.projection = EPSG2WKT(epsg_code)
@property @property
def box(self): def box(self):
mapPoly = get_footprint_polygon(get_corner_coordinates(gt=self.geotransform, cols=self.cols,rows=self.rows)) mapPoly = get_footprint_polygon(get_corner_coordinates(gt=self.geotransform, cols=self.cols,rows=self.rows))
...@@ -237,11 +251,11 @@ class GeoArray(object): ...@@ -237,11 +251,11 @@ class GeoArray(object):
def __getattr__(self, attr): def __getattr__(self, attr):
# check if the requested attribute can not be present because GeoArray has been instanced with an array # check if the requested attribute can not be present because GeoArray has been instanced with an array
if attr not in self.__dict__ and not self.is_inmem and attr in ['shape','dtype','geotransform', 'projection']: if attr not in self.__dir__() and not self.is_inmem and attr in ['shape','dtype','geotransform', 'projection']:
self.set_gdalDataset_meta() self.set_gdalDataset_meta()
if attr in self.__dict__: if attr in self.__dir__(): #__dir__() includes also methods and properties
return self.__dict__[attr] return self.__getattribute__(attr) #__getattribute__ avoids infinite loop
else: else:
raise AttributeError("%s object has no attribute '%s'." %(self.__class__.__name__, attr)) raise AttributeError("%s object has no attribute '%s'." %(self.__class__.__name__, attr))
...@@ -402,24 +416,179 @@ class GeoArray(object): ...@@ -402,24 +416,179 @@ class GeoArray(object):
dill.dump(self,outF) dill.dump(self,outF)
def show(self, xlim=None, ylim=None, band=0, figsize=None, interpolation='none', cmap=None): def _get_plottable_image(self, xlim=None, ylim=None, band=0, res_factor=None, nodataVal=None, out_prj=None):
# handle limits # handle limits
cS, cE = xlim if isinstance(xlim,(tuple,list)) else (0,self.cols-1) cS, cE = xlim if isinstance(xlim, (tuple, list)) else (0, self.cols - 1)
rS, rE = ylim if isinstance(ylim,(tuple,list)) else (0,self.rows-1) rS, rE = ylim if isinstance(ylim, (tuple, list)) else (0, self.rows - 1)
image2plot = self[rS:rE, cS:cE, band]
gt, prj = self.geotransform, self.projection
if res_factor != 1. and image2plot.shape[0] * image2plot.shape[1] > 1e6: # shape > 1000*1000
# sample image down
xdim, ydim = (self.cols * res_factor, self.rows * res_factor) if res_factor else \
tuple((np.array([self.cols, self.rows]) / (np.array([self.cols, self.rows]).max() / 1000))) # normalize
xdim, ydim = int(xdim), int(ydim)
transOpt = ['SRC_METHOD=NO_GEOTRANSFORM'] if tuple(gt)==(0,1,0,0,0,-1) else None
from ...geo.raster.reproject import warp_ndarray
image2plot, gt, prj = warp_ndarray(image2plot, self.geotransform, self.projection,
out_XYdims=(xdim, ydim), in_nodata=nodataVal, out_nodata=nodataVal,
transformerOptions=transOpt, out_prj=out_prj, q=True)
if transOpt and 'NO_GEOTRANSFORM' in ','.join(transOpt):
image2plot = np.flipud(image2plot)
gt=list(gt)
gt[3]=0
print('Note: array has been downsampled to %s x %s for faster visualization.' % (xdim, ydim))
return image2plot, gt, prj
def show(self, xlim=None, ylim=None, band=0, figsize=None, interpolation='none', cmap=None, nodataVal=None,
res_factor=None):
"""Plots the desired array position into a figure.
:param xlim:
:param ylim:
:param band:
:param figsize:
:param interpolation:
:param cmap:
:param nodataVal:
:param res_factor:
:return:
"""
# get image to plot
image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, res_factor, nodataVal)
# set color palette
palette = cmap if cmap else plt.cm.gray
if nodataVal is not None: # do not show nodata
image2plot = np.ma.masked_equal(image2plot, nodataVal)
palette.set_bad('aqua', 0)
palette.set_over ('1')
palette.set_under('0')
# show image # show image
plt.figure(figsize=figsize) plt.figure(figsize=figsize)
plt.imshow(self[rS:rE,cS:cE,band], interpolation=interpolation, cmap=cmap) # FIXME plt.imshow(image2plot, palette, interpolation=interpolation,extent=(0,self.cols,self.rows,0),
#plt.imshow(self[:30,:30,band] if is_3D else self[:30,:30],interpolation='none') vmin=np.percentile(image2plot,2),vmax=np.percentile(image2plot,98),)
plt.show() plt.show()
def show_map(self): def show_map(self, xlim=None, ylim=None, band=0, figsize=None, interpolation='none', cmap=None, nodataVal=None,
raise NotImplementedError res_factor=None, return_map=False):
from mpl_toolkits.basemap import Basemap """
m = Basemap(llcrnrlon=-119, llcrnrlat=22, urcrnrlon=-64, urcrnrlat=49, :param xlim:
projection='lcc', lat_1=33, lat_2=45, lon_0=-95) :param ylim:
:param band:
:param figsize:
:param interpolation:
:param cmap:
:param nodataVal:
:param res_factor:
:param return_map:
:return:
"""
assert self.geotransform and tuple(self.geotransform) != (0,1,0,0,0,-1),\
'A valid geotransform is needed for a map visualization. Got %s.' %self.geotransform
assert self.projection, 'A projection is needed for a map visualization. Got %s.' %self.projection
# get image to plot
image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, res_factor, nodataVal, 'epsg:4326')
# calculate corner coordinates of plot
UL_XY, UR_XY, LR_XY, LL_XY = [(YX[1],YX[0]) for YX in GeoArray(image2plot, gt, prj).box.boxMapYX]
center_lon, center_lat = (UL_XY[0]+UR_XY[0])/2., (UL_XY[1]+LL_XY[1])/2.
# create map
fig = plt.figure(figsize=figsize)
plt.subplots_adjust(left=0.05, right=0.95, top=0.90, bottom=0.05, wspace=0.15, hspace=0.05)
ax = plt.subplot(111)
m = Basemap(projection='tmerc', resolution=None, lon_0=center_lon, lat_0=center_lat,
urcrnrlon=UR_XY[0], urcrnrlat=UR_XY[1], llcrnrlon=LL_XY[0], llcrnrlat=LL_XY[1])
# set color palette
palette = cmap if cmap else plt.cm.gray
if nodataVal is not None: # do not show nodata
image2plot = np.ma.masked_equal(image2plot, nodataVal)
palette.set_bad('aqua', 0)
palette.set_over ('1')
palette.set_under('0')
# add image to map (y-axis must be inverted for basemap)
m.imshow(np.flipud(image2plot), palette, interpolation=interpolation,
vmin=np.percentile(image2plot,2),vmax=np.percentile(image2plot,98))
# add coordinate grid lines
parallels = np.arange(-90, 90., 0.25)
m.drawparallels(parallels, labels=[1, 0, 0, 0], fontsize=12, linewidth=0.4)
meridians = np.arange(-180., 180., 0.25)
m.drawmeridians(meridians, labels=[0, 0, 0, 1], fontsize=12, linewidth=0.4)
if return_map:
return fig,ax, m
else:
plt.show()
def show_map_utm(self, xlim=None, ylim=None, band=0, figsize=None, interpolation='none', cmap=None, nodataVal=None,
res_factor=None, return_map=False):
warnings.warn(UserWarning('This function is still under construction and may not work as expected!'))
# TODO debug this function
# get image to plot
image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, res_factor, nodataVal)
# calculate corner coordinates of plot
box2plot = GeoArray(image2plot, gt, prj).box
UL_XY, UR_XY, LR_XY, LL_XY = [(YX[1], YX[0]) for YX in GeoArray(image2plot, gt, prj).box.boxMapYX]
# Xarr, Yarr = self.box.get_coordArray_MapXY(prj=EPSG2WKT(4326))
UL_XY, UR_XY, LR_XY, LL_XY = [transform_any_prj(self.projection,'epsg:4326',x,y) for y,x in box2plot.boxMapYX]
center_X, center_Y = (UL_XY[0] + UR_XY[0]) / 2., (UL_XY[1] + LL_XY[1]) / 2.
center_lon, center_lat = transform_any_prj(prj,'epsg:4326', center_X, center_Y)
print(center_lon, center_lat)
# create map
fig = plt.figure(figsize=figsize)
plt.subplots_adjust(left=0.05, right=0.95, top=0.90, bottom=0.05, wspace=0.15, hspace=0.05)
ax = plt.subplot(111)
print(UL_XY, UR_XY, LR_XY, LL_XY)
# m = Basemap(projection='tmerc', resolution=None, lon_0=center_lon, lat_0=center_lat,
# urcrnrx=UR_XY[0], urcrnry=UR_XY[1], llcrnrx=LL_XY[0], llcrnry=LL_XY[1])
m = Basemap(projection='tmerc', resolution=None, lon_0=center_lon, lat_