1

I have a single shapefile with many polygons, 100 for this example. I want to iterate through each polygon to generate a Euclidean distance raster for each polygon, using the extents of the initial shapefile. Since each Euclidean distance raster would be run for each individual polygon, isolated, this means that each of these rasters would show increasing distance values spreading out from the borders of the polygon out to the extents, with no other polygons to touch. This means that each of the Euclidean distance rasters would be calculated such that the shapefile only contained a single polygon. I would like to send all of these 100 raster files to the same folder.

How can this be scripted so that individual Euclidean distance raster files are generated for each polygon, where these polygons are treated as being isolated? I am not sure which spatial/raster analysis package would be best and most straightforward for generating these Euclidean distance rasters, whether using python, R, GDAL, etc. Specifically, I am trying to use a for loop so I can perform raster calculations on each produced Euclidean distance raster file, where for this example, I want to add 2 to each raster pixel, such as the following:

for i in polygons_shapefile:
    Euclidean_Distance = generate Euclidean distance(i)
    Altered_Euclidean_Distance = Euclidean_Distance + 2
    Send Altered_Euclidean_Distance to "Outputs" folder

I unfortunately do not have an example shapefile to provide, but really my question could apply to any polygons shapefile.

Update:

I am trying to use the code suggestion from @user2856 with my NYC Boroughs shapefile. For reference, this is what the attribute table of this shapefile looks like: enter image description here

And here is the code I am using to produce a proximity raster of each borough and then store the test_value value:

#Create euclidean distance for each polygon and store "Values"
out_raster_template = "Boroughs_Test/out_{}.tif"
out_proximity_template = "Boroughs_Test/prox_{}.tif"
shape_file = "Boroughs_Test/Boroughs.shp"

pixel_size = 10
nodata = -9999

id_field = 'boro_code'
value_field = 'test_value'

drv = gdal.GetDriverByName("ESRI Shapefile")

shp_ds = gdal.OpenEx(shape_file, gdal.OF_VECTOR)
lyr = shp_ds.GetLayer()

xmin, xmax, ymin, ymax = lyr.GetExtent()
srs = lyr.GetSpatialRef()

feat_def = lyr.GetLayerDefn()

lyr.ResetReading()
for feat in lyr:
    id = feat.GetField(id_field)
    val = feat.GetField(value_field)

    tmp_feat = feat.Clone()

    out_raster = out_raster_template.format(id)
    prox_raster = out_proximity_template.format(id)
    tmp_fn = '/vsimem/tmp.shp'
    tmp_ds = drv.Create(tmp_fn, 0, 0, 0, gdal.GDT_Unknown )
    tmp_lyr = tmp_ds.CreateLayer(tmp_fn, None, feat_def.GetGeomType())
    tmp_lyr.CreateFeature(tmp_feat)
    tmp_feat, tmp_lyr, tmp_ds = None, None, None

    out_ds = gdal.Rasterize(out_raster, tmp_fn,
                   outputType=gdal.GDT_Float32, format='GTIFF', creationOptions=["COMPRESS=DEFLATE"],
                   noData=nodata, initValues=nodata,
                   xRes=pixel_size, yRes=-pixel_size, outputBounds=(xmin, ymin, xmax, ymax), outputSRS=srs,
                   allTouched=True, burnValues=val)

    out_ds = None

    gdal.Translate(prox_raster, out_raster, creationOptions=["COMPRESS=DEFLATE"])
    src_ds = gdal.OpenEx(out_raster, gdal.OF_RASTER)
    dst_ds = gdal.OpenEx(prox_raster, gdal.OF_UPDATE)

    src_band = src_ds.GetRasterBand(1)
    dst_band = dst_ds.GetRasterBand(1)

    gdal.ComputeProximity(src_band, dst_band, options=[f'VALUES={val}'])

    dst_band, src_band, dst_ds, src_ds = None, None, None, None

    drv.Delete(tmp_fn)

This code works and produces the desired rasters, but for some reason prox_1.0.tif and prox_5.0.tif yield no data to show (e.g. a scale of -1.79769e+308 to 1.79769e+308) when viewed in QGIS.

Improve this question
1
  • 2
    A gdal based python solution would look something like open vector dataset, loop through features, gdal.Rasterize each feature, gdal.ComputeProximity each rasterized feature.
    – user2856
    Nov 11, 2021 at 23:50

1 Answer 1

Reset to default
2

Here is an example that

  1. opens vector dataset and loops through each feature,
  2. gdal.Rasterize each feature
  3. gdal.ComputeProximity each rasterized feature

You should reproject your vector if they are in geographic lat/lon CRS so the proximity (distance) values are in metres, I used QGIS to reproject this data to UTM Zone 18N.

from osgeo import gdal
gdal.UseExceptions()

out_raster_template = "out_{}.tif"
shape_file = "boroughs_utm.shp"

pixel_size = 10
nodata = -9999

id_field = 'boro_code'

drv = gdal.GetDriverByName("ESRI Shapefile")

shp_ds = gdal.OpenEx(shape_file, gdal.OF_VECTOR)
lyr = shp_ds.GetLayer()

xmin, xmax, ymin, ymax = lyr.GetExtent()
srs = lyr.GetSpatialRef()

feat_def = lyr.GetLayerDefn()

lyr.ResetReading()
for feat in lyr:
    id = int(feat.GetField(id_field))

    tmp_feat = feat.Clone()

    out_raster = out_raster_template.format(id)
    tmp_fn = '/vsimem/tmp.shp'
    tmp_ds = drv.Create(tmp_fn, 0, 0, 0, gdal.GDT_Unknown )
    tmp_lyr = tmp_ds.CreateLayer(tmp_fn, None, feat_def.GetGeomType())
    tmp_lyr.CreateFeature(tmp_feat)
    tmp_feat, tmp_lyr, tmp_ds = None, None, None

    out_ds = gdal.Rasterize(out_raster, tmp_fn,
                   outputType=gdal.GDT_Float32, format='GTIFF', creationOptions=["COMPRESS=DEFLATE"],
                   noData=nodata, initValues=nodata,
                   xRes=pixel_size, yRes=-pixel_size, outputBounds=(xmin, ymin, xmax, ymax), outputSRS=srs,
                   allTouched=True, burnValues=0)

    out_ds = None

    src_ds = gdal.OpenEx(out_raster, gdal.OF_RASTER)
    dst_ds = gdal.OpenEx(out_raster, gdal.OF_UPDATE)

    src_band = src_ds.GetRasterBand(1)
    dst_band = dst_ds.GetRasterBand(1)

    gdal.ComputeProximity(src_band, dst_band, options=[f'VALUES=0'])

    dst_band, src_band, dst_ds, src_ds = None, None, None, None

    drv.Delete(tmp_fn)

If you want to rasterize some value from your shapefile to use later, but still do the euc distance to another:

from osgeo import gdal
gdal.UseExceptions()


#Create euclidean distance for each polygon and store "Values"
out_raster_template = "out_{}.tif"
out_proximity_template = "prox_{}.tif"
shape_file = "boroughs_utm.shp"

pixel_size = 10
nodata = -9999

id_field = 'boro_code'
value_field = 'value'

drv = gdal.GetDriverByName("ESRI Shapefile")

shp_ds = gdal.OpenEx(shape_file, gdal.OF_VECTOR)
lyr = shp_ds.GetLayer()

xmin, xmax, ymin, ymax = lyr.GetExtent()
srs = lyr.GetSpatialRef()

feat_def = lyr.GetLayerDefn()

lyr.ResetReading()
for feat in lyr:
    id = int(feat.GetField(id_field))
    val = feat.GetField(value_field)

    tmp_feat = feat.Clone()

    out_raster = out_raster_template.format(id)
    prox_raster = out_proximity_template.format(id)
    tmp_fn = '/vsimem/tmp.shp'
    tmp_raster = '/vsimem/tmp.tif'
    tmp_ds = drv.Create(tmp_fn, 0, 0, 0, gdal.GDT_Unknown )
    tmp_lyr = tmp_ds.CreateLayer(tmp_fn, None, feat_def.GetGeomType())
    tmp_lyr.CreateFeature(tmp_feat)
    tmp_feat, tmp_lyr, tmp_ds = None, None, None

    out_ds = gdal.Rasterize(out_raster, tmp_fn,
                            outputType=gdal.GDT_Float32, format='GTIFF', creationOptions=["COMPRESS=DEFLATE"],
                            noData=nodata, initValues=nodata,
                            xRes=pixel_size, yRes=-pixel_size, outputBounds=(xmin, ymin, xmax, ymax), outputSRS=srs,
                            allTouched=True, burnValues=val)

    out_ds = None

    out_ds = gdal.Rasterize(tmp_raster, tmp_fn,
                            outputType=gdal.GDT_Int32, format='GTIFF', creationOptions=["COMPRESS=DEFLATE"],
                            noData=nodata, initValues=nodata,
                            xRes=pixel_size, yRes=-pixel_size, outputBounds=(xmin, ymin, xmax, ymax), outputSRS=srs,
                            allTouched=True, burnValues=id)

    out_ds = None

    gdal.Translate(prox_raster, out_raster, creationOptions=["COMPRESS=DEFLATE"])
    src_ds = gdal.OpenEx(tmp_raster, gdal.OF_RASTER)
    dst_ds = gdal.OpenEx(prox_raster, gdal.OF_UPDATE)

    src_band = src_ds.GetRasterBand(1)
    dst_band = dst_ds.GetRasterBand(1)

    gdal.ComputeProximity(src_band, dst_band, options=[f'VALUES={id}'])

    dst_band, src_band, dst_ds, src_ds = None, None, None, None

    drv.Delete(tmp_fn)

enter image description here

Improve this answer
10
  • Thanks for the suggestion! Why are you using a temp folder instead of just a usual directory folder? I reprojected the shapefile as you suggested, but would I then need to send it to a folder called "tmp"? I tried just using the original directory for the input shapefile, but got the error: AttributeError: 'NoneType' object has no attribute 'GetRasterBand', pointing to an error with this line: ---> 44 src_band = src_ds.GetRasterBand(1). Also, does out_raster_template also need to be in this "tmp" folder?
    – LostinSpatialAnalysis
    Nov 28, 2021 at 2:26
  • Put your data anywhere you like. No temp folder in above code (now). Either use a full path to your data or set your working directory appropriately and use a relative path.
    – user2856
    Nov 28, 2021 at 3:05
  • 'NoneType' object has no attribute 'GetRasterBand' means GDAL couldn't open your raster and returned None instead of a gdal.Dataset instead of raising an exception, usually because the path is wrong and you then tried to GetRasterBand on None. If you add gdal.UseExceptions() to your script (after importing gdal), then you will get a more useful error.
    – user2856
    Nov 28, 2021 at 3:44
  • Thanks for these explanations, I was able to successfully implement your first code suggestion above, where I created a separate Euclidean distance raster for each borough polygon. However, for some reason, I am having an issue with the second code suggestion, where I create a separate Euclidean distance raster for each borough polygon, and store a different value as another rasterized layer to use later, which you labelled as 'VALUE'. Though, running this produces out_1.0, out_2.0, out_3.0, and out_5.0, missing 4.0, and prox_1.0, prox_2.0, and prox_3.0, missing 4.0 and 5.0.
    – LostinSpatialAnalysis
    Nov 30, 2021 at 1:53
  • 1
    I just ran the code copied from your Q and the issue was that the values are floating point, I adjusted code so it uses the ID field to compute proximity to.
    – user2856
    Jan 26 at 23:36

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.