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?

  • What code snippet did you test? What was the precise result of running it? – PolyGeo Dec 8 '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
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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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