I have a geo-referenced tiff file and a point shapefile. I want to ignore all the points in the shapefile which are outside raster's boundary.
Asked
Active
Viewed 284 times
0
-
Duplicate? https://gis.stackexchange.com/questions/57964/get-vector-features-inside-a-specific-extent – RoperMaps Jun 21 '17 at 13:05
-
@GeoWill the answer to that post is bit unclear to me, how did he manage to calculate the bounding box, there is not much explanation. – FJ_Abbasi Jun 21 '17 at 13:36
-
Either scroll down to the second answer which references this method of getting raster extent. Or open your raster in qgis and look at the properties. – RoperMaps Jun 21 '17 at 14:17
-
what is your meaning of "ignore" all the points in the shapefile which are outside raster's boundary? – xunilk Jun 22 '17 at 02:28
-
@xunilk deleting the points. – FJ_Abbasi Jun 22 '17 at 08:37
-
Based in your more recent commentary, I think that is preferable to produce a new point vector layer, with points overlay raster, for preserving original point layer. Please, see my answer. – xunilk Jun 22 '17 at 12:34
1 Answers
1
Based in your more recent commentary, I think that is preferable to produce a new point vector layer, with points overlay raster, for preserving original point layer. Next PyQGIS code produces it as memory layer in QGIS. You only have to change my name layers for yours.
registry = QgsMapLayerRegistry.instance()
raster = registry.mapLayersByName('noise_linear_raster')
vector = registry.mapLayersByName('random_points3')
extent_polygon = raster[0].extent().asWktPolygon()
geom_extent = QgsGeometry.fromWkt(extent_polygon)
feat_vector = [ feat for feat in vector[0].getFeatures() ]
new_feat_vector = [ feat for feat in feat_vector if feat.geometry().within(geom_extent) ]
epsg = vector[0].crs().postgisSrid()
uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"
mem_layer = QgsVectorLayer(uri,
'no_ignored_points',
'memory')
prov = mem_layer.dataProvider()
for i, feat in enumerate(new_feat_vector):
feat.setAttributes([i])
prov.addFeatures(new_feat_vector)
QgsMapLayerRegistry.instance().addMapLayer(mem_layer)
I tried out above code in Python Console of QGIS and I got:
Red points are corresponding to memory layer and all of them are within raster layer.
xunilk
- 29,891
- 4
- 41
- 80