3

I have some tifs files inside a folder and want to check if a shape coordinates are inside one of the tifs with pyqgis.

I tried this from Determining coordinates of corners of raster layer using PyQGIS? to get the raster corners so I can check if my shape is inside it :

from osgeo import gdal, osr

bag = gdal.Open('LC82220762017084LGN00_B2.TIF')  # replace it with your file # raster is projected
bag_gtrn = bag.GetGeoTransform()
bag_proj = bag.GetProjectionRef()
bag_srs = osr.SpatialReference(bag_proj)
geo_srs =bag_srs.CloneGeogCS()                 # new srs obj to go 
from x,y -> φ,λ
transform = osr.CoordinateTransformation( bag_srs, geo_srs)

bag_bbox_cells = (
(0., 0.),
(0, bag.RasterYSize),
(bag.RasterXSize, bag.RasterYSize),
(bag.RasterXSize, 0),
)

geo_pts = []
for x, y in bag_bbox_cells:
    x2 = bag_gtrn[0] + bag_gtrn[1] * x + bag_gtrn[2] * y
    y2 = bag_gtrn[3] + bag_gtrn[4] * x + bag_gtrn[5] * y
    geo_pt = transform.TransformPoint(x2, y2)[:2]
    geo_pts.append(geo_pt)
    print x, y, '->', x2, y2, '->', geo_pt

But it returned

SyntaxError: Non-ASCII character '\xcf' in file intersect.py on line 8, but no encoding declared; see http://python.org/dev/peps/pep-0263/ for details

And tried this Determining if shapefile and raster overlap in Python using OGR/GDAL? to check if my shape overlaps the raster:

import ogr, gdal

raster = gdal.Open('LC82220762017084LGN00_B2.TIF')
vector = ogr.Open('Faz_Esperanca.shp')

# Get raster geometry
transform = raster.GetGeoTransform()
pixelWidth = transform[1]
pixelHeight = transform[5]
cols = raster.RasterXSize
rows = raster.RasterYSize

xLeft = transform[0]
yTop = transform[3]
xRight = xLeft+cols*pixelWidth
yBottom = yTop-rows*pixelHeight

ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(xLeft, yTop)
ring.AddPoint(xLeft, yBottom)
ring.AddPoint(xRight, yTop)
ring.AddPoint(xRight, yBottom)
ring.AddPoint(xLeft, yTop)
rasterGeometry = ogr.Geometry(ogr.wkbPolygon)
rasterGeometry.AddGeometry(ring)

layer = vector.GetLayer()
feature = layer.GetFeature(0)
vectorGeometry = feature.GetGeometryRef()

print rasterGeometry.Intersect(vectorGeometry)

Running this with my rasters and vector I get False returning every time even if the vector overlaps the raster or not.

Are there any function on Qgis that I not used?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Frox
  • 31
  • 6
  • What code snippet did you test? What was the precise result of running it? – PolyGeo Dec 08 '17 at 20:23
  • Thanks for the help. I tried to edit the question with some results of running the accepted answers of the questions that I checked before posting. – Frox Dec 11 '17 at 14:33
  • Welcome to GIS SE! As a new user be sure to take the [Tour]. For questions that involve code we ask that you show us where you are stuck with your own code by including a code snippet in your question. There is an [edit] button beneath your question which will enable you to do that and a {} button that enables you to format any highlighted code nicely. We need to know precisely where you are stuck with one test that you run and include a code snippet for within your question body. – PolyGeo Dec 11 '17 at 22:45
  • Thanks for the help, I edited it again and put the codes that I tried to run. – Frox Dec 12 '17 at 10:54
  • Which of the two tests do you wish to ask about in this particular question? – PolyGeo Dec 12 '17 at 11:10
  • In fact, none of them, my question is if there is a way to check if vector and raster overlaps in PyQgis – Frox Dec 12 '17 at 11:35
  • The only answer on either of those questions that seems to use PyQGIS is https://gis.stackexchange.com/a/57824/115 so you could try focusing your question on that. – PolyGeo Dec 12 '17 at 18:59

1 Answers1

0

After days of researching I resolved this using Postgis. 

SELECT ST_Intersects( ST_GeographyFromText('SRID=4326;POLYGON((-51.829893 -22.063562,-50.015800 -22.435339,-50.417151 -24.163293,-52.254290 -23.786699,-51.829893 -22.063562))'),ST_GeographyFromText('SRID=4326;POINT(-51.22224 -23.67402)'))

It returns true if POINT is inside POLYGON.

I use a loop to check if my shape borders are inside tif borders.

This the only code that I find which works with EPSG:4326 (Latitude Longitude, spheroidal values)

Frox
  • 31
  • 6