conversion.py 4.83 KB
Newer Older
1
2
3
4
5
# -*- coding: utf-8 -*-

import numpy as np

# custom
Daniel Scheffler's avatar
Daniel Scheffler committed
6
7
from shapely.geometry import shape, mapping, box
from shapely.geometry import Polygon  # noqa F401  # flake8 issue
8

Daniel Scheffler's avatar
Daniel Scheffler committed
9
10
11
12
13
14
15
16
try:
    from osgeo import ogr
    from osgeo import osr
    from osgeo import gdal
except ImportError:
    import ogr
    import osr
    import gdal
17
18

from ..coord_trafo import imYX2mapYX, mapYX2imYX, pixelToMapYX
19
from ...dtypes.conversion import get_dtypeStr, dTypeDic_NumPy2GDAL
Daniel Scheffler's avatar
Daniel Scheffler committed
20

21
22
__author__ = "Daniel Scheffler"

23

24
25
26
27
28
29
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])
30
31
32
33
34
35
36
37
38
39
40
    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)


41
def shapelyBox2BoxYX(shapelyBox, coord_type='image'):
42
    xmin, ymin, xmax, ymax = shapelyBox.bounds
43
44
45
    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))
46
47
48
49
    return UL_YX, UR_YX, LR_YX, LL_YX


def get_boxImXY_from_shapelyPoly(shapelyPoly, im_gt):
Daniel Scheffler's avatar
Daniel Scheffler committed
50
    # type: (Polygon, tuple) -> list
51
52
53
54
55
    """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
    """
56
57
58
59
    def get_coordsArr(shpPoly): return 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]
60
61
62
    return boxImXY


63
def round_shapelyPoly_coords(shapelyPoly, precision=10, out_dtype=None):
64
65
66
67
    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
68
69
70
    return shape(geojson)


Daniel Scheffler's avatar
Daniel Scheffler committed
71
72
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
73
74
75
76
77
78
79
    """
    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
80
    :param fillVal: fill value used to fill in where no point geometry is available
Daniel Scheffler's avatar
Daniel Scheffler committed
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
118
    """

    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
119
120
    if fillVal is not None:
        band.Fill(fillVal)
Daniel Scheffler's avatar
Daniel Scheffler committed
121
122
123
124
125
126
    band.FlushCache()

    # Rasterize
    gdal.RasterizeLayer(target_ds, [1], layer, options=["ATTRIBUTE=VAL"])

    out_arr = target_ds.GetRasterBand(1).ReadAsArray()
127
    out_gt = target_ds.GetGeoTransform()
Daniel Scheffler's avatar
Daniel Scheffler committed
128
129
    out_prj = target_ds.GetProjection()

130
    del target_ds, ds, layer, band
Daniel Scheffler's avatar
Daniel Scheffler committed
131
132

    return out_arr, out_gt, out_prj