conversion.py 4.82 KB
Newer Older
1
2
3
4
5
6
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"

import numpy as np

# custom
Daniel Scheffler's avatar
Daniel Scheffler committed
7
from shapely.geometry import shape, mapping, box, Polygon
Daniel Scheffler's avatar
Daniel Scheffler committed
8
9
10
11
12
13
14
15
try:
    from osgeo import ogr
    from osgeo import osr
    from osgeo import gdal
except ImportError:
    import ogr
    import osr
    import gdal
16
17
18

from ..coord_trafo import imYX2mapYX, mapYX2imYX, pixelToMapYX

Daniel Scheffler's avatar
Daniel Scheffler committed
19
20
from py_tools_ds.ptds.dtypes.conversion import get_dtypeStr, dTypeDic_NumPy2GDAL

21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48


def shapelyImPoly_to_shapelyMapPoly_withPRJ(shapelyImPoly, gt,prj):
    #ACTUALLY PRJ IS NOT NEEDED BUT THIS FUNCTION RETURNS OTHER VALUES THAN shapelyImPoly_to_shapelyMapPoly
    geojson   = mapping(shapelyImPoly)
    coords    = list(geojson['coordinates'][0])
    coordsYX = pixelToMapYX(coords,geotransform=gt,projection=prj)
    coordsXY = tuple([(i[1],i[0]) for i in coordsYX])
    geojson['coordinates'] = (coordsXY,)
    return shape(geojson)


def shapelyImPoly_to_shapelyMapPoly(shapelyBox, gt):
    xmin, ymin, xmax, ymax = shapelyBox.bounds
    ymax, xmin = imYX2mapYX((ymax, xmin), gt)
    ymin, xmax = imYX2mapYX((ymin, xmax), gt)
    return box(xmin, ymin, xmax, ymax)


def shapelyBox2BoxYX(shapelyBox,coord_type='image'):
    xmin, ymin, xmax, ymax = shapelyBox.bounds
    assert coord_type in ['image','map']
    UL_YX, UR_YX, LR_YX, LL_YX = ((ymin,xmin),(ymin,xmax),(ymax,xmax),(ymax,xmin)) if coord_type=='image' else \
                                 ((ymax,xmin),(ymax,xmax),(ymin,xmax),(ymin,xmin))
    return UL_YX, UR_YX, LR_YX, LL_YX


def get_boxImXY_from_shapelyPoly(shapelyPoly, im_gt):
Daniel Scheffler's avatar
Daniel Scheffler committed
49
    # type: (Polygon, tuple) -> list
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
    """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.
    :param shapelyPoly:     <shapely.Polygon>
    :param im_gt:           <list> the GDAL geotransform of the target image
    """
    get_coordsArr = lambda shpPoly: np.swapaxes(np.array(shpPoly.exterior.coords.xy), 0, 1)
    coordsArr     = get_coordsArr(shapelyPoly)
    boxImXY       = [mapYX2imYX((Y, X), im_gt) for X, Y in coordsArr.tolist()] # FIXME incompatible to GMS version
    boxImXY       = [(i[1],i[0]) for i in boxImXY]
    return boxImXY


def round_shapelyPoly_coords(shapelyPoly, precision=10,out_dtype=None):
    geojson = mapping(shapelyPoly)
    geojson['coordinates'] = np.round(np.array(geojson['coordinates']), precision)
    if out_dtype:
        geojson['coordinates'] = geojson['coordinates'].astype(out_dtype)
Daniel Scheffler's avatar
Daniel Scheffler committed
67
68
69
    return shape(geojson)


Daniel Scheffler's avatar
Daniel Scheffler committed
70
71
def points_to_raster(points, values, tgt_res, prj=None, fillVal=None):
    # type: (np.ndarray, np.ndarray, float, str, float) -> (np.ndarray, list, str)
Daniel Scheffler's avatar
Daniel Scheffler committed
72
73
74
75
76
77
78
    """
    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
Daniel Scheffler's avatar
Daniel Scheffler committed
79
    :param fillVal: fill value used to fill in where no point geometry is available
Daniel Scheffler's avatar
Daniel Scheffler committed
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
    """

    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)
Daniel Scheffler's avatar
Daniel Scheffler committed
118
119
    if fillVal is not None:
        band.Fill(fillVal)
Daniel Scheffler's avatar
Daniel Scheffler committed
120
121
122
123
124
125
126
127
128
129
130
131
132
133
    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