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

added functions

dtypes.conversion:
- added get_dtypeStr() (moved here from CoReg_Sat)

geo.vector.conversion:dict_to_proj4dict_to_proj4
- added points_to_raster(): new rasterization function

geo.projection:
- added dict_to_proj4()

- updated __version__
parent 209ca2f2
......@@ -15,7 +15,7 @@ __all__=[#'compatibility',
'similarity',
'GeoArray']
__version__ = '20170306_01'
__version__ = '20170315_01'
__author__='Daniel Scheffler'
# Validate GDAL version
......
......@@ -2,6 +2,8 @@
__author__ = "Daniel Scheffler"
import datetime
import numpy as np
try:
from osgeo import gdal
......@@ -43,3 +45,11 @@ dTypeDic_NumPy2GDALcompatible = \
dict(zip(dTypeDic_NumPy2GDAL.keys(),
[dTypeDic_GDAL2Numpy[dTypeDic_NumPy2GDAL[str(np.dtype(NDT))]] for NDT in dTypeDic_NumPy2GDAL.keys()]))
def get_dtypeStr(val):
is_numpy = 'numpy' in str(type(val))
DType = str(np.dtype(val)) if is_numpy else 'int' if isinstance(val,int) else 'float' if isinstance(val,float) else \
'str' if isinstance(val,str) else 'complex' if isinstance(val,complex) else \
'date' if isinstance(val,datetime.datetime) else None
assert DType, 'data type not understood'
return DType
\ No newline at end of file
......@@ -4,6 +4,7 @@ __author__ = "Daniel Scheffler"
import re
import os
import pyproj
# custom
try:
......@@ -50,6 +51,14 @@ def proj4_to_dict(proj4):
return {k:v for k,v in [kv.split('=') for kv in items]}
def dict_to_proj4(proj4dict):
# type: (dict) -> str
"""Converts a PROJ4-like dictionary into a PROJ4 string.
:param proj4dict: <dict> the PROJ4-like dictionary
"""
return pyproj.Proj(proj4dict).srs
def prj_equal(prj1, prj2):
#type: (any,any) -> bool
"""Checks if the given two projections are equal.
......
......@@ -5,9 +5,19 @@ import numpy as np
# custom
from shapely.geometry import shape, mapping, box
try:
from osgeo import ogr
from osgeo import osr
from osgeo import gdal
except ImportError:
import ogr
import osr
import gdal
from ..coord_trafo import imYX2mapYX, mapYX2imYX, pixelToMapYX
from py_tools_ds.ptds.dtypes.conversion import get_dtypeStr, dTypeDic_NumPy2GDAL
def shapelyImPoly_to_shapelyMapPoly_withPRJ(shapelyImPoly, gt,prj):
......@@ -55,3 +65,66 @@ def round_shapelyPoly_coords(shapelyPoly, precision=10,out_dtype=None):
if out_dtype:
geojson['coordinates'] = geojson['coordinates'].astype(out_dtype)
return shape(geojson)
def points_to_raster(points, values, tgt_res, prj=None):
# type: (np.ndarray, np.ndarray, float, str) -> np.ndarray, list, str
"""
Converts a set of point geometries with associated values into a raster array.
:param points: list or 1D numpy.ndarray containings shapely.geometry point geometries
:param values: list or 1D numpy.ndarray containing int or float values
:param tgt_res: target resolution in projection units
:param prj: WKT projection string
"""
ds = ogr.GetDriverByName("Memory").CreateDataSource('wrk')
if prj is not None:
srs = osr.SpatialReference()
srs.ImportFromWkt(prj)
else:
srs = None
layer = ds.CreateLayer('', srs, ogr.wkbPoint)
# create field
DTypeStr = get_dtypeStr(values if isinstance(values, np.ndarray) else values[0])
FieldType = ogr.OFTInteger if DTypeStr.startswith('int') else ogr.OFTReal
FieldDefn = ogr.FieldDefn('VAL', FieldType)
if DTypeStr.startswith('float'):
FieldDefn.SetPrecision(6)
layer.CreateField(FieldDefn) # Add one attribute
for i in range(len(points)):
# Create a new feature (attribute and geometry)
feat = ogr.Feature(layer.GetLayerDefn())
feat.SetGeometry(ogr.CreateGeometryFromWkb(points[i].wkb)) # Make a geometry, from Shapely object
feat.SetField('VAL', values[i])
layer.CreateFeature(feat)
feat.Destroy()
x_min, x_max, y_min, y_max = layer.GetExtent()
# Create the destination data source
cols = int((x_max - x_min) / tgt_res)
rows = int((y_max - y_min) / tgt_res)
target_ds = gdal.GetDriverByName('MEM').Create('raster', cols, rows, 1, dTypeDic_NumPy2GDAL[DTypeStr])
target_ds.SetGeoTransform((x_min, tgt_res, 0, y_max, 0, -tgt_res))
target_ds.SetProjection(prj if prj else '')
band = target_ds.GetRasterBand(1)
band.FlushCache()
# Rasterize
gdal.RasterizeLayer(target_ds, [1], layer, options=["ATTRIBUTE=VAL"])
out_arr = target_ds.GetRasterBand(1).ReadAsArray()
out_gt = target_ds.GetGeoTransform()
out_prj = target_ds.GetProjection()
target_ds = ds = layer = band = None
return out_arr, out_gt, out_prj
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