diff --git a/src/core/accessibility/__init__.py b/src/core/accessibility/__init__.py new file mode 100644 index 0000000..40a96af --- /dev/null +++ b/src/core/accessibility/__init__.py @@ -0,0 +1 @@ +# -*- coding: utf-8 -*- diff --git a/src/core/accessibility/create_grid.py b/src/core/accessibility/create_grid.py new file mode 100644 index 0000000..d056c0c --- /dev/null +++ b/src/core/accessibility/create_grid.py @@ -0,0 +1,100 @@ +# -*- coding: utf-8 -*- +""" +/*************************************************************************** + + GeoHealth + A QGIS plugin + + ------------------- + begin : 2016-11-20 + copyright : (C) 2014 by Etienne Trimaille + email : etienne@trimaille.eu + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ +""" + +from PyQt4.QtCore import QVariant + +from qgis.core import ( + QGis, + QgsVectorLayer, + QgsCoordinateReferenceSystem, + QgsPoint, + QgsFields, + QgsField, + QgsFeature, + QgsGeometry, +) + +from GeoHealth.src.core.tools import create_memory_layer + + +def create_grid(extent, space, crs): + """Create a grid polygon. + + :param extent: The extend where to create the grid. + :type extent: QgsRectangle + + :param space: The X and Y size. + :type space: float + + :param crs: The CRS of the new grid. + :type crs: QgsCoordinateReferenceSystem + + :return: THe vector grid. + :rtype: QgsVectorLayer + """ + + fields = QgsFields() + fields.append(QgsField('id', QVariant.Int, '', 10, 0)) + fields.append(QgsField('xmin', QVariant.Double, '', 24, 15)) + fields.append(QgsField('xmax', QVariant.Double, '', 24, 15)) + fields.append(QgsField('ymin', QVariant.Double, '', 24, 15)) + fields.append(QgsField('ymax', QVariant.Double, '', 24, 15)) + + writer = create_memory_layer('grid', QGis.Polygon, crs, fields) + writer.startEditing() + + feat = QgsFeature() + feat.initAttributes(len(fields)) + feat.setFields(fields) + geom = QgsGeometry() + id_polygon = 0 + + count = 0 + y = extent.yMaximum() + while y >= extent.yMinimum(): + x = extent.xMinimum() + while x <= extent.xMaximum(): + polygon = [ + [ + QgsPoint(x, y), + QgsPoint(x + space, y), + QgsPoint(x + space, y - space), + QgsPoint(x, y - space), + QgsPoint(x, y) + ] + ] + feat.setGeometry(geom.fromPolygon(polygon)) + feat.setAttribute(0, id_polygon) + feat.setAttribute(1, x) + feat.setAttribute(2, x + space) + feat.setAttribute(3, y - space) + feat.setAttribute(4, y) + writer.addFeature(feat) + id_polygon += 1 + x += space + y -= space + count += 1 + + writer.commitChanges() + + return writer diff --git a/src/core/accessibility/process.py b/src/core/accessibility/process.py new file mode 100644 index 0000000..2ba9749 --- /dev/null +++ b/src/core/accessibility/process.py @@ -0,0 +1,213 @@ +def create_grid(size): + """Create a polygonal grid using Processing. + + :param size: The cell size. + :type size: int + + :return: The grid layer in memory. + :rtype: QgsVectorLayer + """ + output_filename = unique_filename(prefix='grid', suffix='.shp') + + result = processing.runalg( + 'qgis:vectorgrid', + '336199.970553,352338.397991,7636164.67975,7648562.41208', + size, # X spacing + size, # Y spacing + 0, # Output as polygons + output_filename) + + layer = QgsVectorLayer(output_filename, 'grid', 'ogr') + layer.setCrs(QgsCoordinateReferenceSystem(32740)) + + remove_fields(layer, ['xmin', 'xmax', 'ymin', 'ymax']) + + # Make a copy in memory + memory = create_memory_layer( + 'grid', layer.geometryType(), layer.crs(), layer.fields()) + copy_layer(layer, memory) + + print "NB cells : %s" % layer.featureCount() + + return memory + + +def intersect_grid_and_roads(grid, roads): + """Create a subset of the grid.""" + + grid.startEditing() + grid.addAttribute( + QgsField('has_road', QVariant.String, len=10, prec=0)) + intersection_index = grid.fieldNameIndex('has_road') + + spatial_index = create_spatial_index(roads) + + for cell in grid.getFeatures(): + geometry = cell.geometry() + intersects = spatial_index.intersects(geometry.boundingBox()) + + for i in intersects: + request = QgsFeatureRequest().setFilterFid(i) + road = next(roads.getFeatures(request)) + road_geometry = road.geometry() + + if geometry.intersects(road_geometry): + has_intersection = True + break + else: + has_intersection = False + + grid.changeAttributeValue( + cell.id(), intersection_index, has_intersection) + + grid.commitChanges() + return grid + + +def centroids(layer): + """Create centroids.""" + centroids_layer = create_memory_layer( + 'centroids', + QGis.Point, + layer.crs(), + layer.fields()) + + centroids_layer.startEditing() + + for feature in layer.getFeatures(): + new_feature = QgsFeature(feature) + new_feature.setGeometry(feature.geometry().centroid()) + centroids_layer.addFeature(new_feature) + centroids_layer.commitChanges() + return centroids_layer + + +def create_graph(network, centroid_layer, target_point_layer): + """Create a graph based on a road network and some tied points.""" + tied_points = [] + for f in centroid_layer.getFeatures(): + tied_points.append(f.geometry().centroid().asPoint()) + for feature in target_point_layer.getFeatures(): + tied_points.append(feature.geometry().asPoint()) + + graph = Graph(network, tied_points) + return graph + + +def intersecting_blocks(line, destination, grid): + """Function to fetch intersectings polygons from the grid.""" + dest_id_field = grid.fieldNameIndex('destination_id') + request = QgsFeatureRequest() + request.setFilterRect(line.boundingBox()) + request.setFilterExpression('"destination_id" is None') + for feature in grid.getFeatures(request): + if feature.geometry().intersects(line): + grid.changeAttributeValue(feature.id(), dest_id_field, destination) + + +def assign_cost_to_cells(network_graph, source, destination, id_field): + """Assign the nearest destination point from the source layer. + + :param network_graph: The network graph. + :type network_graph: Graph + + :param source: Grid as a polygon vector layer + :type source: QgsVectorLayer + + :param destination: The destination point layer. + :type destination: QgsVectorLayer + + :param id_field: The ID field in the destination layer. + :type id_field: basestring + """ + spatial_index = create_spatial_index(destination) + destination_features = {} + for feature in destination.getFeatures(): + destination_features[feature.id()] = feature + index_id_field = destination.fieldNameIndex(id_field) + + fields = [QgsField('distance', QVariant.Int)] + routes = create_memory_layer( + 'routes', QGis.Line, destination.crs(), fields) + routes.startEditing() + + source.startEditing() + source.addAttribute( + QgsField('destination_id', QVariant.Int, len=5, prec=0)) + dest_id_field = source.fieldNameIndex('destination_id') + source.addAttribute(QgsField('distance', QVariant.Int, len=5, prec=0)) + distance_field = source.fieldNameIndex('distance') + source.commitChanges() + + request = QgsFeatureRequest() + request.setFilterExpression('"has_road" = \'true\'') + # request.setLimit(20) # Hack for now to speedup development + i = 0 + for source_cell in source.getFeatures(request): + source_geometry_point = source_cell.geometry().centroid().asPoint() + desination_id = source_cell['destination_id'] + if desination_id is None or isinstance(desination_id, QPyNullVariant): + source.startEditing() + nearest_health_points = spatial_index.nearestNeighbor( + source_geometry_point, 5) + minimum_distance = None + minimal_geom = None + for health_point in nearest_health_points: + try: + i += 1 + p = destination_features[health_point].geometry().asPoint() + geom, distance, _ = network_graph.route( + source_geometry_point, p) + except: + distance = -1 + + if minimum_distance is None or ( + minimum_distance > distance >= 0): + minimal_geom = geom + minimum_distance = distance + + if minimum_distance: + feature = QgsFeature() + feature.setGeometry(minimal_geom) + feature.setAttributes([minimum_distance]) + routes.addFeatures([feature]) + + index_id = destination_features[health_point][index_id_field] + destination_value = index_id + else: + destination_value = '-1' + minimum_distance = '-1' + + source.changeAttributeValue( + source_cell.id(), dest_id_field, destination_value) + source.changeAttributeValue( + source_cell.id(), distance_field, minimum_distance) + intersecting_blocks(minimal_geom, destination_value, grid) + source.commitChanges() + + else: + print 'speedup' + geom, distance, _ = network_graph.route( + source_geometry_point, + destination_features[desination_id].geometry().asPoint()) + + source.commitChanges() + # source.commitErrors() + routes.commitChanges() + print 'Call : %s' % i + return routes, source + +if health_points.crs != roads.crs(): + health_points = reproject(health_points, roads.crs()) + +grid = create_grid(cell_size) +# show_qgis_layer(grid) +grid = intersect_grid_and_roads(grid, roads) +# show_qgis_layer(cells_without_roads) +# centroids_layer = centroids(cells_with_roads) +# show_qgis_layer(centroids_layer) +network = create_graph(roads, grid, health_points) +routes, grids = assign_cost_to_cells(network, grid, health_points, id_field) +show_qgis_layer(grids) +show_qgis_layer(routes) +show_qgis_layer(health_points) \ No newline at end of file diff --git a/src/core/accessibility/test/test_grid.py b/src/core/accessibility/test/test_grid.py new file mode 100644 index 0000000..8361018 --- /dev/null +++ b/src/core/accessibility/test/test_grid.py @@ -0,0 +1,32 @@ +import unittest + +from qgis.core import QgsRectangle, QgsCoordinateReferenceSystem +from qgis.testing.mocked import get_iface +from qgis.utils import iface + +from GeoHealth.src.core.accessibility.create_grid import create_grid + +if iface: + APP = iface +else: + APP = get_iface() + + +class TestTest(unittest.TestCase): + + def setUp(self): + pass + + def tearDown(self): + pass + + def test_grid(self): + + extent = QgsRectangle(0, 0, 4, 4) + size = 1 + crs = QgsCoordinateReferenceSystem(3857) + print crs.authid() + self.assertTrue(crs.isValid()) + layer = create_grid(extent, 2, crs) + print layer.featureCount() + print "TOTO" diff --git a/src/doc/help.py b/src/doc/help.py index 5e44873..c080d4e 100644 --- a/src/doc/help.py +++ b/src/doc/help.py @@ -326,6 +326,23 @@ def help_incidence_point(): return html +def help_accessibility(): + title = tr('Accessibility') + intro = tr('To write this later.') + inputs = [ + tr('To write inputs 1'), + tr('To write inputs 2')] + outputs = [ + tr('Output 1'), + tr('Output 2'), + ] + more = [ + tr('More information how we do that'), + ] + html = html_table(title, intro, inputs, outputs, more) + return html + + def help_attribute_table(): title = tr('Export attribute table') intro = tr('Export as CSV format without geometry.') diff --git a/src/gui/analysis/accessibility.py b/src/gui/analysis/accessibility.py new file mode 100644 index 0000000..1b06722 --- /dev/null +++ b/src/gui/analysis/accessibility.py @@ -0,0 +1,40 @@ +# -*- coding: utf-8 -*- +""" +/*************************************************************************** + + GeoHealth + A QGIS plugin + + ------------------- + begin : 2014-08-20 + copyright : (C) 2014 by Etienne Trimaille + email : etienne@trimaille.eu + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ +""" + +from PyQt4.QtGui import ( + QWidget, QDialogButtonBox, QFileDialog, QApplication) +from PyQt4.QtCore import pyqtSignal, QSettings, QVariant + +from GeoHealth.src.utilities.resources import get_ui_class + +FORM_CLASS = get_ui_class('analysis', 'accessibility.ui') + + +class AccessibilityWidget(QWidget, FORM_CLASS): + + signalAskCloseWindow = pyqtSignal(name='signalAskCloseWindow') + + def __init__(self, parent=None): + self.parent = parent + super(AccessibilityWidget, self).__init__() + self.setupUi(self) diff --git a/src/gui/main_window.py b/src/gui/main_window.py index e690806..000fa24 100644 --- a/src/gui/main_window.py +++ b/src/gui/main_window.py @@ -35,6 +35,7 @@ help_incidence, help_incidence_point, help_density_point, + help_accessibility, help_attribute_table, ) from GeoHealth.src.gui.import_gui.open_shapefile import ( @@ -43,6 +44,7 @@ from GeoHealth.src.gui.import_gui.open_raster import OpenRasterWidget from GeoHealth.src.gui.import_gui.open_xls_dbf import ( OpenXlsDbfFileWidget) +from GeoHealth.src.gui.analysis.accessibility import AccessibilityWidget from GeoHealth.src.gui.analysis.blur_dialog import BlurWidget from GeoHealth.src.gui.analysis.stats_dialog import StatsWidget from GeoHealth.src.gui.analysis.incidence_dialog import IncidenceDialog @@ -177,6 +179,13 @@ def __init__(self, parent=None): } } ] + }, { + 'label': 'Accessibility', + 'icon': resource('accessibility.svg'), + 'content': { + 'widget': AccessibilityWidget(), + 'help': help_accessibility() + } } ] }, diff --git a/src/resources/accessibility.svg b/src/resources/accessibility.svg new file mode 100644 index 0000000..13c4ad2 --- /dev/null +++ b/src/resources/accessibility.svg @@ -0,0 +1,21 @@ + + + + + + + + \ No newline at end of file diff --git a/src/ui/analysis/accessibility.ui b/src/ui/analysis/accessibility.ui new file mode 100644 index 0000000..e804de1 --- /dev/null +++ b/src/ui/analysis/accessibility.ui @@ -0,0 +1,297 @@ + + + Accessibility + + + + 0 + 0 + 651 + 695 + + + + Incidence + + + + + + QFormLayout::ExpandingFieldsGrow + + + + + Network + + + + + + + + 0 + 0 + + + + + + + + Extent + + + + + + + + 0 + 0 + + + + + + + + Cell size + + + + + + + + 0 + 0 + + + + 99 + + + + + + + Points + + + + + + + + 0 + 0 + + + + + + + + NB points + + + + + + + + 0 + 0 + + + + + + + + Output + + + + + + + + + Save to temporary directory + + + + + + + Browse + + + + + + + + + + + + + + Add a symbology + + + true + + + false + + + + + + + + Farest + + + + + + + + 255 + 246 + 246 + + + + + + + + Qt::Horizontal + + + + 40 + 20 + + + + + + + + Closest + + + + + + + + 202 + 33 + 36 + + + + + + + + + + + + Classes + + + + + + + 1 + + + 5 + + + + + + + Qt::Horizontal + + + + 40 + 20 + + + + + + + + Mode + + + + + + + + + + + + + + + Qt::Horizontal + + + QDialogButtonBox::Cancel|QDialogButtonBox::Ok + + + + + + + Qt::Vertical + + + + 20 + 40 + + + + + + + + + QgsColorButtonV2 + QToolButton +
qgis.gui
+
+ + QgsMapLayerComboBox + QComboBox +
qgis.gui
+
+ + QgsCollapsibleGroupBox + QGroupBox +
qgis.gui
+ 1 +
+
+ + +