Commit 8d78e630 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

some improvements

geo.vector.conversion:
- points_to_raster(): added keyword fillVal

geo.vector.geometry:
- boxObj.boundsIm: coordinates are now rounded to integers

geo.projection:
- added proj4_to_WKT()

io.raster.GeoArray:
- GeoArray.show_map(): added keyword 'ax'
- _clip_array_at_mapPos(): bugfix

- updated __version__
parent e2a126ff
...@@ -15,7 +15,7 @@ __all__=[#'compatibility', ...@@ -15,7 +15,7 @@ __all__=[#'compatibility',
'similarity', 'similarity',
'GeoArray'] 'GeoArray']
__version__ = '20170316_01' __version__ = '20170328_01'
__author__='Daniel Scheffler' __author__='Daniel Scheffler'
# Validate GDAL version # Validate GDAL version
......
...@@ -59,6 +59,16 @@ def dict_to_proj4(proj4dict): ...@@ -59,6 +59,16 @@ def dict_to_proj4(proj4dict):
return pyproj.Proj(proj4dict).srs return pyproj.Proj(proj4dict).srs
def proj4_to_WKT(proj4str):
# type: (str) -> str
"""Converts a PROJ4-like string into a WKT string.
:param proj4str: <dict> the PROJ4-like string
"""
srs = osr.SpatialReference()
srs.ImportFromProj4(proj4str)
return srs.ExportToWkt()
def prj_equal(prj1, prj2): def prj_equal(prj1, prj2):
#type: (any,any) -> bool #type: (any,any) -> bool
"""Checks if the given two projections are equal. """Checks if the given two projections are equal.
......
...@@ -4,7 +4,7 @@ __author__ = "Daniel Scheffler" ...@@ -4,7 +4,7 @@ __author__ = "Daniel Scheffler"
import numpy as np import numpy as np
# custom # custom
from shapely.geometry import shape, mapping, box from shapely.geometry import shape, mapping, box, Polygon
try: try:
from osgeo import ogr from osgeo import ogr
from osgeo import osr from osgeo import osr
...@@ -46,7 +46,7 @@ def shapelyBox2BoxYX(shapelyBox,coord_type='image'): ...@@ -46,7 +46,7 @@ def shapelyBox2BoxYX(shapelyBox,coord_type='image'):
def get_boxImXY_from_shapelyPoly(shapelyPoly, im_gt): def get_boxImXY_from_shapelyPoly(shapelyPoly, im_gt):
# type: (shapely.Polygon,tuple) -> np.ndarray # type: (Polygon, tuple) -> list
"""Converts each vertex coordinate of a shapely polygon into image coordinates corresponding to the given """Converts each vertex coordinate of a shapely polygon into image coordinates corresponding to the given
geotransform without respect to invalid image coordinates. Those must be filtered later. geotransform without respect to invalid image coordinates. Those must be filtered later.
:param shapelyPoly: <shapely.Polygon> :param shapelyPoly: <shapely.Polygon>
...@@ -67,8 +67,8 @@ def round_shapelyPoly_coords(shapelyPoly, precision=10,out_dtype=None): ...@@ -67,8 +67,8 @@ def round_shapelyPoly_coords(shapelyPoly, precision=10,out_dtype=None):
return shape(geojson) return shape(geojson)
def points_to_raster(points, values, tgt_res, prj=None): def points_to_raster(points, values, tgt_res, prj=None, fillVal=None):
# type: (np.ndarray, np.ndarray, float, str) -> np.ndarray, list, str # type: (np.ndarray, np.ndarray, float, str, float) -> (np.ndarray, list, str)
""" """
Converts a set of point geometries with associated values into a raster array. Converts a set of point geometries with associated values into a raster array.
...@@ -76,6 +76,7 @@ def points_to_raster(points, values, tgt_res, prj=None): ...@@ -76,6 +76,7 @@ def points_to_raster(points, values, tgt_res, prj=None):
:param values: list or 1D numpy.ndarray containing int or float values :param values: list or 1D numpy.ndarray containing int or float values
:param tgt_res: target resolution in projection units :param tgt_res: target resolution in projection units
:param prj: WKT projection string :param prj: WKT projection string
:param fillVal: fill value used to fill in where no point geometry is available
""" """
ds = ogr.GetDriverByName("Memory").CreateDataSource('wrk') ds = ogr.GetDriverByName("Memory").CreateDataSource('wrk')
...@@ -114,6 +115,8 @@ def points_to_raster(points, values, tgt_res, prj=None): ...@@ -114,6 +115,8 @@ def points_to_raster(points, values, tgt_res, prj=None):
target_ds.SetGeoTransform((x_min, tgt_res, 0, y_max, 0, -tgt_res)) target_ds.SetGeoTransform((x_min, tgt_res, 0, y_max, 0, -tgt_res))
target_ds.SetProjection(prj if prj else '') target_ds.SetProjection(prj if prj else '')
band = target_ds.GetRasterBand(1) band = target_ds.GetRasterBand(1)
if fillVal is not None:
band.Fill(fillVal)
band.FlushCache() band.FlushCache()
# Rasterize # Rasterize
......
...@@ -14,7 +14,7 @@ from ..projection import prj_equal ...@@ -14,7 +14,7 @@ from ..projection import prj_equal
class boxObj(object): class boxObj(object):
def __init__(self, **kwargs): def __init__(self,**kwargs):
"""Create a dynamic/self-updating box object that represents a rectangular or quadratic coordinate box """Create a dynamic/self-updating box object that represents a rectangular or quadratic coordinate box
according to the given keyword arguments. according to the given keyword arguments.
Note: Either mapPoly+gt or imPoly+gt or wp+ws or boxMapYX+gt or boxImYX+gt must be passed. Note: Either mapPoly+gt or imPoly+gt or wp+ws or boxMapYX+gt or boxImYX+gt must be passed.
...@@ -112,7 +112,7 @@ class boxObj(object): ...@@ -112,7 +112,7 @@ class boxObj(object):
@property @property
def boundsIm(self): def boundsIm(self):
"""Returns xmin,xmax,ymin,ymax in image coordinates.""" """Returns xmin,xmax,ymin,ymax in image coordinates."""
boxImXY = get_boxImXY_from_shapelyPoly(self.mapPoly, self.gt) boxImXY = [(round(x, 0), round(y, 0)) for x, y in get_boxImXY_from_shapelyPoly(self.mapPoly, self.gt)]
return corner_coord_to_minmax(boxImXY) # xmin,xmax,ymin,ymax return corner_coord_to_minmax(boxImXY) # xmin,xmax,ymin,ymax
@property @property
......
...@@ -950,15 +950,16 @@ class GeoArray(object): ...@@ -950,15 +950,16 @@ class GeoArray(object):
plt.show() plt.show()
def show_map(self, xlim=None, ylim=None, band=0, boundsMap=None, boundsMapPrj=None, figsize=None, def show_map(self, xlim=None, ylim=None, band=0, boundsMap=None, boundsMapPrj=None, ax=None, figsize=None,
interpolation='none', cmap=None, nodataVal=None, res_factor=None, return_map=False, zoomable=False): interpolation='none', cmap=None, nodataVal=None, res_factor=None, return_map=False, zoomable=False):
""" """
:param xlim: :param xlim:
:param ylim: :param ylim:
:param band: :param band: band index (starting with 0)
:param boundsMap: xmin, ymin, xmax, ymax :param boundsMap: xmin, ymin, xmax, ymax
:param boundsMapPrj: :param boundsMapPrj:
:param ax: allows to pass a matplotlib axis object where figure is plotted into
:param figsize: :param figsize:
:param interpolation: :param interpolation:
:param cmap: :param cmap:
...@@ -1004,7 +1005,7 @@ class GeoArray(object): ...@@ -1004,7 +1005,7 @@ class GeoArray(object):
# create map # create map
fig = plt.figure(figsize=figsize) 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) 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) ax = ax if ax is not None else plt.subplot(111)
m = Basemap(projection='tmerc', resolution=None, lon_0=center_lon, lat_0=center_lat, 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]) urcrnrlon=UR_XY[0], urcrnrlat=UR_XY[1], llcrnrlon=LL_XY[0], llcrnrlat=LL_XY[1])
...@@ -1021,16 +1022,20 @@ class GeoArray(object): ...@@ -1021,16 +1022,20 @@ class GeoArray(object):
palette.set_under('0') palette.set_under('0')
# add image to map (y-axis must be inverted for basemap) # add image to map (y-axis must be inverted for basemap)
#vmin=None # TODO make this adjustable
#vmax=None
if zoomable: if zoomable:
m.imshow(image2plot, palette, interpolation=interpolation, vmin=vmin, vmax=vmax) m.imshow(image2plot, palette, interpolation=interpolation, vmin=vmin, vmax=vmax)
else: else:
m.imshow(np.flipud(image2plot), palette, interpolation=interpolation, vmin=vmin, vmax=vmax) m.imshow(np.flipud(image2plot), palette, interpolation=interpolation, vmin=vmin, vmax=vmax)
# add coordinate grid lines # add coordinate grid lines
parallels = np.arange(-90, 90., 0.25) parallels = np.arange(-90, 90., 0.25) # TODO make this adjustable
#parallels = np.arange(-90, 90., 0.1)
m.drawparallels(parallels, labels=[1, 0, 0, 0], fontsize=12, linewidth=0.4) m.drawparallels(parallels, labels=[1, 0, 0, 0], fontsize=12, linewidth=0.4)
meridians = np.arange(-180., 180., 0.25) meridians = np.arange(-180., 180., 0.25)
#meridians = np.arange(-180., 180., 0.1)
m.drawmeridians(meridians, labels=[0, 0, 0, 1], fontsize=12, linewidth=0.4) m.drawmeridians(meridians, labels=[0, 0, 0, 1], fontsize=12, linewidth=0.4)
if return_map: if return_map:
...@@ -1458,10 +1463,10 @@ def _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=None, fillVal=0): ...@@ -1458,10 +1463,10 @@ def _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=None, fillVal=0):
out_gt[0], out_gt[3] = xmin, ymax out_gt[0], out_gt[3] = xmin, ymax
# get image area to read # get image area to read
cS, rS = [int(i) for i in mapXY2imXY((xmin, ymax), arr_gt)] cS, rS = [int(i) for i in mapXY2imXY((xmin, ymax), arr_gt)] # UL
cE, rE = [int(i)-1 for i in mapXY2imXY((xmax, ymin), arr_gt)] cE, rE = [int(i)-1 for i in mapXY2imXY((xmax, ymin), arr_gt)] # LR
if 0 <= rS <= rows - 1 and 0 <= rE <= rows - 1 and 0 <= cS <= cols - 1 and 0 <= cE <= rows - 1: if 0 <= rS <= rows - 1 and 0 <= rE <= rows - 1 and 0 <= cS <= cols - 1 and 0 <= cE <= cols - 1:
"""requested area is within the input array""" """requested area is within the input array"""
if bands==1: if bands==1:
out_arr = arr[rS:rE + 1, cS:cE + 1] out_arr = arr[rS:rE + 1, cS:cE + 1]
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment