Python >> Tutoriel Python >  >> Python Tag >> Array

Lire un fichier de formes sous forme de tableau en utilisant Python ?

Votre script était généralement correct, mais vous n'avez pas modifié le nom du champ attributaire que vous vouliez pixelliser.

Dans l'exemple que vous avez posté, vous définissez ['ATTRIBUTE=ID'] en tant que champ, mais il n'existe pas dans votre fichier de formes. Vous n'avez que "Habitats" et "surface" en tant que champs, vous devez donc modifier correctement le code.

Par conséquent, vous deviez modifier les dossiers des calques de fichiers de formes et rastérisés, ainsi que le crs.

J'ai légèrement modifié le code de cette manière :

import gdal
from osgeo import osr
from osgeo import ogr

def layer(shapefile):

    # 1) opening the shapefile
    source_ds = ogr.Open(shapefile)
    source_layer = source_ds.GetLayer()

    # 2) Creating the destination raster data source

    pixelWidth = pixelHeight = 1 # depending how fine you want your raster
    x_min, x_max, y_min, y_max = source_layer.GetExtent()
    cols = int((x_max - x_min) / pixelHeight)
    rows = int((y_max - y_min) / pixelWidth)
    target_ds = gdal.GetDriverByName('GTiff').Create(raster_path, cols, rows, 1, gdal.GDT_Byte) 
    target_ds.SetGeoTransform((x_min, pixelWidth, 0, y_min, 0, pixelHeight))
    band = target_ds.GetRasterBand(1)
    NoData_value = 255
    band.SetNoDataValue(NoData_value)
    band.FlushCache()

    # 4) Instead of setting a general burn_value, use optionsand set it to the attribute that contains the relevant unique value ["ATTRIBUTE=ID"]
    gdal.RasterizeLayer(target_ds, [1], source_layer, options = ['ATTRIBUTE=surface'])

    # 5) Adding a spatial reference
    target_dsSRS = osr.SpatialReference()
    target_dsSRS.ImportFromEPSG(2975)
    target_ds.SetProjection(target_dsSRS.ExportToWkt())
    return gdal.Open(raster_path).ReadAsArray()


raster_path = 'C:/Users/path_to_the_rasterized_output/temp.tif'

shapefile = 'C:/Users/path_to_the_shapefile/shapefile_maido_tipe.shp'

print layer(shapefile)

et je pense que cela fonctionne maintenant car j'obtiens ce calque pixellisé (qui chevauche le fichier de formes):

et ce retour du print layer(shapefile) ligne (vous ne voyez que la valeur '255' car vous l'avez définie comme valeur nodata):

[[255 255 255 ..., 255 255 255]
 [255 255 255 ..., 255 255 255]
 [255 255 255 ..., 255 255 255]
 ..., 
 [255 255 255 ..., 255 255 255]
 [255 255 255 ..., 255 255 255]
 [255 255 255 ..., 255 255 255]]