0

Earlier I asked a question, about Python GeoPandas crops result or sets boundary to data

And I still cant understand how to work with geospatial data.

I generate coordinates and create a plot with this code:

import csv

import numpy
from matplotlib import pyplot
from matplotlib.tri import Triangulation


def read(path):
    result = list()

    with open(path) as csv_handler:
        csv_reader = csv.reader(csv_handler)

        for line in csv_reader:
            for item in line[0].split(';'):
                item = float(item)

                if item == -999.0:
                    result.append(numpy.nan)
                elif item == -999.7999877929688:
                    result.append(numpy.nan)
                else:
                    result.append(item)

    return result


def nan_equal(actual):
    try:
        numpy.testing.assert_equal(actual, numpy.nan)
    except AssertionError:
        return False

    return True


def clear(longitudes, latitudes, values):
    if len(longitudes) == len(latitudes) == len(values):
        result_longitudes = list()
        result_latitudes = list()
        result_values = list()

        for i in range(len(longitudes)):
            if not nan_equal(longitudes[i]):
                if not nan_equal(latitudes[i]):
                    if not nan_equal(values[i]):
                        result_longitudes.append(longitudes[i])
                        result_latitudes.append(latitudes[i])
                        result_values.append(values[i])

        return result_longitudes, result_latitudes, result_values
    else:
        raise ValueError


def main():
    size = 5000

    longitudes = numpy.random.rand(size)
    latitudes = numpy.random.rand(size)
    values = numpy.random.uniform(-1000.0, 1000.0, size)

    triangulation1 = Triangulation(longitudes, latitudes)
    triangulation2 = Triangulation(longitudes, latitudes)
    triangulation3 = Triangulation(longitudes, latitudes)

    def set_mask(triangulation, alpha=0.4):
        triangles = triangulation.triangles

        xtri = longitudes[triangles] - numpy.roll(longitudes[triangles], 1, axis=1)
        ytri = latitudes[triangles] - numpy.roll(latitudes[triangles], 1, axis=1)
        maxi = numpy.max(numpy.sqrt(xtri ** 2 + ytri ** 2), axis=1)

        triangulation.set_mask(maxi > alpha)

    set_mask(triangulation2, alpha=0.1)
    set_mask(triangulation3, alpha=0.3)

    figure, (axes1, axes2, axes3) = pyplot.subplots(ncols=3, figsize=(60, 20))

    axes1.tricontourf(triangulation1, values)
    axes1.scatter(longitudes, latitudes, s=5, color='black')

    axes2.tricontourf(triangulation2, values)
    axes2.scatter(longitudes, latitudes, s=5, color='black')

    axes3.tricontourf(triangulation3, values)
    axes3.scatter(longitudes, latitudes, s=5, color='black')

    pyplot.show()


if __name__ == '__main__':
    main()

But if i change random coordinates generation to:

longitudes = numpy.random.uniform(-180.0, 180.0, size)
latitudes = numpy.random.uniform(-90.0, 90.0, size)

Code doesnt work and write error message

ValueError: z array must not contain non-finite values within the triangulation

Update

I found some code to create polygon by points, but I don't know how to convert it to triangular mask

def alpha_shape(points, alpha):
    if len(points) < 4:
        return geometry.MultiPoint(list(points)).convex_hull

    def add_edge(edges, edge_points, coords, i, j):
        if (i, j) in edges or (j, i) in edges:
            return
        edges.add((i, j))
        edge_points.append(coords[[i, j]])

    coords = numpy.array([point.coords[0] for point in points])

    tri = Delaunay(coords)
    edges = set()
    edge_points = []
    for ia, ib, ic in tri.vertices:
        pa = coords[ia]
        pb = coords[ib]
        pc = coords[ic]

        a = math.sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
        b = math.sqrt((pb[0] - pc[0]) ** 2 + (pb[1] - pc[1]) ** 2)
        c = math.sqrt((pc[0] - pa[0]) ** 2 + (pc[1] - pa[1]) ** 2)

        s = (a + b + c) / 2.0

        area = math.sqrt(s * (s - a) * (s - b) * (s - c))
        circum_r = a * b * c / (4.0 * area)

        if circum_r < 1.0 / alpha:
            add_edge(edges, edge_points, coords, ia, ib)
            add_edge(edges, edge_points, coords, ib, ic)
            add_edge(edges, edge_points, coords, ic, ia)

    m = geometry.MultiLineString(edge_points)
    triangles = list(polygonize(m))

    concave_hull = cascaded_union(triangles)

    return concave_hull, edge_points 
3
  • 2
    this seems like an XY problem - could you describe what you are trying to do as well as what you have tried so far.
    – Ian Turton
    Apr 14, 2021 at 7:37
  • 2
    As I told you in your previous question, It is a matplotlib problem, not a geospatial problem.
    – gene
    Apr 14, 2021 at 7:59
  • 2
    "I find some code to create polygon by points": where ? (concave_hulls.ipynb for example).
    – gene
    Apr 14, 2021 at 8:28

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.