Tag Archives: python

Using GDAL to get raster extent

Hi there guys!!!

Let’s suppose we want to determine the extent of a raster file and we want to use GDAL and Python. How can we do that?

Let’s start importing GDAL:

from osgeo import gdal

Then we need t0 open the raster file:

gdalSrc = gdal.Open('/foo/bar/filename')

To finish getting what we need, let’s get our affine transform coefficients with the following:

upx, xres, xskew, upy, yskew, yres = gdalSrc.GetGeotransform()

Where upx is the upper pixel x, xres is the pixel width and xskew is the shear in the x direction as we can see in the following picture (https://en.wikipedia.org/wiki/Affine_transformation):

xskew

Let’s go forward. The variable upy is the upper pixel y, yres is the pixel height and yskew is the shear in the y direction as we can see in the following picture (https://en.wikipedia.org/wiki/Affine_transformation):

yskew

Understanting the variables we can use the documentation obtained at http://www.gdal.org/gdal_datamodel.html to get the following relationship:

Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)

Where GT(2) is the xskew and GT(4) is the yskew.

With all of this together, we can make the following code snippet to accomplish our mission:

from osgeo import gdal, ogr
gdalSrc = gdal.Open('/foo/bar/filename')
upx, xres, xskew, upy, yskew, yres = gdalSrc.GetGeotransform()
cols = gdalSrc.RasterXSize
rows = gdalSrc.RasterYSize

ulx = upx + 0*xres + 0*xskew
uly = upy + 0*yskew + 0*yres

llx = upx + 0*xres + rows*xskew
lly = upy + 0*yskew + rows*yres

lrx = upx + cols*xres + rows*xskew
lry = upy + cols*yskew + rows*yres

urx = upx + cols*xres + 0*xskew
ury = upy + cols*yskew + 0*yres

I hope this can be useful to you guys.

Advertisements

Deaggregate geometries with pyqgis

Have you ever needed to explode multi geometry layer into a single geometry layer, using in each new geometry the attributes of the original multi one?  If you were working with FME, for instance, you basically would just use the transformer Deaggregator. Let’s learn how to solve this problem with python and QGIS!

The following code snippet teaches you how to work with QgsVectorLayers, it’s attributes and how to manipulate geometries.

from qgis.core import QgsVectorLayer, QgsFeature, QgsMapLayerRegistry

#fill in your input layer name. In this example, our inputLyrName is input_layer
inputLyrName = 'input_layer'
inputLyr = QgsMapLayerRegistry.instance().mapLayersByName(inputLyrName)[0]

#fill in your output layer name. In this example, our outputLyrName is output_layer
outputName = 'output_layer'
outputLyr = QgsMapLayerRegistry.instance().mapLayersByName(outputLyrName)[0]

#tests type of output: if it is a multi parted geometry or
#a single parted geometry
if outputLyr.wkbType() in [QGis.WKBPoint, QGis.WKBLineString, QGis.WKBPolygon]:
     isMulti = False
else:
     isMulti = True

outputLyr.startEditing()
addList = []
for feat in inputLyr.getFeatures():
     #gets all parts of geometry as an individual single geometry
     parts = feat.geometry().asGeometryCollection()
     #checks if it isMulti, if it is, convert each
     #part in geometryCollection to multi
     if isMulti:
          for part in parts:
               part.convertToMultiType()

     #for each part, get original set of attribute and create a new feat
     #with this set
     for i in range(0,len(parts)):
          #new feature constructor. newFeat has all atributes of feat
          newFeat = QgsFeature(feat)
          #set geometry with part
          newFeat.setGeometry(parts[i])
          #get field id and get defaultValue from provider
          idx = newFeat.fieldNameIndex('id')
          newFeat.setAttribute(idx,provider.defaultValue(idx))
          addList.append(newFeat)
outputLyr.addFeatures(addList,True)
outputLyr.commitChanges()

Hope you guys like it!

Using GRASS with pyqgis to clean up geometries

We all know that GRASS is a great GIS software. Combined with QGIS it is even more great!

Using GRASS from within QGIS is very useful to deal with daily GIS problems. Everyone that works with geospatial data knows how annoying is to clean up geometries full of errors. The manual process demands lots of time and we can always forget something in the end. Do this kind of job automatically is faster and safer.

Let’s se how to do this using pyqgis. Imagine that we have a database layer like this:

error

A good way to clean problems like those shown above and at the same time solve snapping problems is to use the following tools in v.clean.advanced provided by GRASS:

  • break
  • rmsa
  • rmdangle

If you want a description on how those tools work, take a look at: https://grass.osgeo.org/grass73/manuals/v.clean.html

To clean it using grass we can use the following piece of code:


#choosing the algorithm
alg = 'grass7:v.clean.advanced'

#getting the vector layer we want to clean
input = iface.activeLayer()

#setting tools
tools = 'break,rmsa,rmdangle'
threshold = -1
#getting mapcanvas extent (bounding box) supposing we can see our data
e = iface.mapCanvas().extent()
xmax = e.xMaximum()
ymax = e.yMaximum()
xmin = e.xMinimum()
ymin = e.yMinimum()

extent = '{0},{1},{2},{3}'.format(xmin, xmax, ymin, ymax)

#setting parameters: choose them according to your data
snap = 100.0
minArea = 0.001

#running the grass algorithm
ret = processing.runalg(alg, input, tools, threshold, extent, snap, minArea, None, None)

#getting output layer
outputLayer = processing.getObject(ret['output'])
#Adding to registry
QgsMapLayerRegistry.instance().addMapLayer(outputLayer)

#getting error flags
errorLayer = processing.getObject(ret['error'])
#Adding to registry
QgsMapLayerRegistry.instance().addMapLayer(errorLayer) 

After running, we can see results like this:

cleaned

Quite good, right?