0

I'm trying to develop a plugin for QGIS where I want to extrapolate points to a surface. The plugin needs to handle 100,000 points rather quickly and following code works alright. The problem is the extent.

file_name&file_name_poly = string with the path to the files.
shpfile = path to shapefile with the points
extent = '%f,%f,%f,%f'% (shpfile.bbox[0],shpfile.bbox[2],
                         shpfile.bbox[1],shpfile.bbox[3])
processing.runalg("grass:v.voronoi", file_name + ".shp", False, True,
                  extent, -1, 0, 3, file_name_poly)

Here is a picture of the result, as you can see, the polygons to the left looks rather alright, but the polygons to right becomes huge. Detailed picture of the points and their polygon

So I wonder is it possible to set a maximum size for each polygon? Or is it somehow possible to set the extent as a polygon?

Edit!

I have tried to implement this solution hence the code became much slower and with the same problem. This picture is a part of the field, the dots are the original data points, the solid polygons are the ones from the grass processing voronoi algorithm and the dotted line polygons are from ST_DelaunayTriangles link.

New voronoi picture This just a part of field, I wasn't able to run the complete field, since it was taking too long time. Hence, the problem remains that I can't set an outer boundary of the field.

Writing this post might solved the problem. The solution might be to do the grass voronoi diagram first and the make an st_within to remove the unnecessary parts afterwards.

I'll post again when I've tried that solution!

1 Answer 1

0

Well, I solved the problem with the st_intersect function, the implementation is not perfect it is a little bit slow hence it works! Here is a part of the code that does the "magic". mplPath is part of matplolib. defined_field is the field as the user manual choose. shp is shapefile imported as shp. The polygon_file is the output from the voronoi processing. temp_polygon is temporary postgres table no_intersecting_point is a function that checks that the point is within "defined_field". This works not 100% because of the projection so I "cheated" and added an ST_buffer on the defined_field.

def no_intersecting_point(self, shape):
    field_shape = mplPath.Path(self.defined_field)
    for point in shape.points:
        if field_shape.contains_point([point[0], point[1]]):
            return False
    return True


with shp.Reader(polgyon_file) as polygon:
    shapes = polygon.shapes()
    sql_raw = "Insert into temp_polygon(polygon) VALUES "
    for i in range(len(shapes)):
        points = ''
        if self.no_intersecting_point(shapes[i]):
            continue
        for j in range(len(shapes[i].points)):
            points += str(shapes[i].points[j][0]) + ' ' + str(shapes[i].points[j][1]) + ", "
        sql_raw += "(ST_Intersection(ST_buffer(ST_GeomFromText('POLYGON(" + str(self.defined_field)[1:-1].replace(',',' ').replace(')  (', ',') +")', 4326),0.0001) ,ST_GeomFromText('POLYGON((" + points[:-2] + "))', 4326))), "
    sql = sql_raw[:-2]
    self.DB.insert_data_polygon(sql, tbl_name)

Hope this can help someone.

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.