GDAL has a pretty handy proximity calculator where each pixel in the output raster represents the distance to the nearest pixel in the input raster. For instance, distance to the nearest road. Is there any transformer doing that in FME ? Thanks!
Best answer by fiorinoh
Hi @fiorinoh,
I think, we can do it using standard FME tools. Here is how I did it - first, I generate a distance raster where each pixel represents the distance between the center pixel and all other pixels (I use a point cloud approach for this). Then, I turn all the road network pixels to points and pass my distance raster to them. Then, I mosaic all rasters together using minimal pixel values for overlapping rasters. As a last step, I deal with pixels beyond Distance parameter setting them to NoData. I tried to do this with the original distance raster, but somehow, it significantly affects performance (we'll investigate).
I used the road network of the city where I live. Once rasterized (1m resolution), it creates over 144,000 pixels representing roads.
My laptop generates the output in 1 minute if I set distance to 200 meters (400x400 distance raster), in 2.5 minutes with 400 meters (800x800 pixels), and under 4 minutes for 500 m distance (1000x1000 pixels). The final raster is about 8K pixels along each side.
I'll see whether we can compare it somehow with the GDAL proximity calculator. I also can make a custom transformer from this workspace if you think this method does what you expect and its performance is reasonable. See the attached template. It was made with FME 2022, but it opens and runs with FME 2021.2 as well.
Dmitri
Hi @dmitribagh , this is my last solution to this issue. As it is not easy to compute proximity rasters with FME transformers, I decided to use a python caller + GDAL. That was not easy because FME raster type and GDAL raster type are not compatible. Hope this code will help the FME community. BR.
import fmeobjects
import numpy as np
from osgeo import gdal
class MyBandTilePopulator(fmeobjects.FMEBandTilePopulator):
"""
This is a subclass of the FMEBandTilePopulator superclass.
It will be used when data is requested to create a new tile and
and populate it to a new FMEBand.
"""
def __init__(self, rasterData):
self._rasterData = rasterData
# required method
def clone(self):
"""
This method is used to create a copy of the data
multiple times while creating a new band
"""
return MyBandTilePopulator(self._rasterData)
# required method
def getTile(self, startRow, startCol, tile):
"""
Creates a new tile that's sized based on the input tile.
Populates that tile using this populator's raster data beginning
at the startRow and startCol.
"""
numRows, numCols = tile.getNumRows(), tile.getNumCols()
newTile = fmeobjects.FMEUInt16Tile(numRows, numCols)
data = newTile.getData()
for row in range(startRow, startRow+numRows):
for col in range(startCol, startCol+numCols):
if row < len(self._rasterData) and col < len(self._rasterData[0]):
data[row - startRow][col - startCol] = self._rasterData[row][col]
newTile.setData(data)
return newTile
class FeatureProcessor(object):
def __init__(self):
self.rasterData = []
def input(self, feature):
self.sysRef = feature.getCoordSys()
raster = feature.getGeometry()
rp = raster.getProperties()
band = raster.getBand(0)
bp = band.getProperties()
print("=== FME Input Raster ===")
print("Coordinate System = " + self.sysRef)
self.pixelWidth = rp.getSpacingX()
print("pixelWidth = " + str(self.pixelWidth))
self.pixelHeight = rp.getSpacingY()
print("pixelHeight = " + str(self.pixelHeight))
self.originX = rp.getOriginX()
print("originX = " + str(self.originX))
self.originY = rp.getOriginY()
print("originY = " + str(self.originY))
numRows = rp.getNumRows()
print("numRows = " + str(numRows))
numCols = rp.getNumCols()
print("numCols = " + str(numCols))
print("===")
tile = fmeobjects.FMEGray16Tile(numRows, numCols)
bandData = band.getTile(0, 0, tile).getData()
array = np.array(bandData)
driver = gdal.GetDriverByName('GTiff')
num_of_bands = 1
print("Creating file raster_transportation_roads.tiff")
srcRaster = driver.Create('raster_transportation_roads.tiff', numCols, numRows, num_of_bands, gdal.GDT_UInt16)
srcRaster.SetGeoTransform((self.originX, self.pixelWidth, 0, self.originY, 0, self.pixelHeight))
srcBand = srcRaster.GetRasterBand(1)
srcBand.WriteArray(array)
srcBand.FlushCache()
srcBand = None
srcRaster = None
src_ds = gdal.Open('raster_transportation_roads.tiff')
geotransform = src_ds.GetGeoTransform()
srcBand = src_ds.GetRasterBand(1)
print("Origin = ({}, {})".format(geotransform[0], geotransform[3]))
print("Pixel Size = ({}, {})".format(geotransform[1], geotransform[5]))
print("===")
print("Creating file proximity_transportation_roads.tiff with GDAL")
dst_filename='proximity_transportation_roads.tiff'
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create(dst_filename, numCols, numRows, num_of_bands, gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform(geotransform)
dst_ds.SetProjection(src_ds.GetProjectionRef())
dstBand = dst_ds.GetRasterBand(1)
gdal.ComputeProximity(srcBand, dstBand, ["DISTUNITS=PIXEL"])
self.rasterData = dstBand.ReadAsArray().tolist()
dstBand.FlushCache()
dstBand = None
srcBand = None
print("===")
def close(self):
# creating the raster data and specifying the formatting of the new raster
rasterData = self.rasterData
# specifying all of the properties for the new FMERaster
numRows, numCols = len(rasterData), len(rasterData[0])
xCellOrigin, yCellOrigin = 0.5, 0.5
xSpacing, ySpacing = self.pixelWidth, self.pixelHeight
xOrigin, yOrigin = self.originX, self.originY
xRotation, yRotation = 0.0, 0.0
# creating the new FMERaster
rasterProperties = fmeobjects.FMERasterProperties(numRows, numCols,
xSpacing, ySpacing,
xCellOrigin, yCellOrigin,
xOrigin, yOrigin,
xRotation, yRotation)
raster = fmeobjects.FMERaster(rasterProperties)
# Populating the contents of the band and appending it to the raster
bandTilePopulator = MyBandTilePopulator(rasterData)
bandName = ''
bandProperties = fmeobjects.FMEBandProperties(bandName,
fmeobjects.FME_INTERPRETATION_UINT16,
fmeobjects.FME_TILE_TYPE_FIXED,
numRows, numCols)
band = fmeobjects.FMEBand(bandTilePopulator, rasterProperties,
bandProperties)
raster.appendBand(band)
# creating a new feature with the FMERaster geometry to be output
feature = fmeobjects.FMEFeature()
feature.setGeometry(raster)
feature.setCoordSys(self.sysRef)
self.pyoutput(feature)
def process_group(self):
"""When 'Group By' attribute(s) are specified, this method is called
once all the FME Features in a current group have been sent to input().
FME Features sent to input() should generally be cached for group-by
processing in this method when knowledge of all Features is required.
The resulting Feature(s) from the group-by processing should be emitted
through self.pyoutput().
FME will continue calling input() a number of times followed
by process_group() for each 'Group By' attribute, so this
implementation should reset any class members for the next group.
"""
pass
+19- Safer
Hello @fiorinoh, good question. I had a look for equivalent tools in FME, but I’m not sure we have a transformer that 100% fits the bill (so-to-speak). While looking into this, I found a script on the GDAL (gdal_proximity.py) website you may be able to use with the PythonCaller.
For an alternative method, you might be able to use point clouds! You could try using a RasterExpressionEvaluator to extract the raster cells of interest. Once those particular cells are isolated, you could convert them to a pointcloud with a PointCloudCombiner. From there you might be able to find a clever way to calculate distance between your points (PointCloudExpressionEvaluator maybe). I think this will take a lot of trial/error though. I was getting ideas from this blog (viewshed analysis section).
There was also a suggestion I saw to use the RasterConvolver, but I think it uses kernels (eg. 3x3) so it would be per pixel, not distance. There is an option to use a custom kernel, but I think you'd need to have a really in-depth understanding of your data to use this.
Sorry I don’t have one clear solution to present you! Hopefully this will give you somewhere to start though! Maybe someone else will have a better idea! Best, Kailin.
Thank you @kailinatsafe . Yes I tried with PythonCaller, it is fine you can import gdal and use the ComputeProximity function etc. The problem is FME Raster type is not compatible with GDAL Raster type (as far as I know...) Therefore, you have to put all the raw data of FME objects into GDAL objects. No impossible but ugh! Above all, it misses the point of a nice and tidy FME workspace: why not use GDAL for the whole geotreatment? Though, I have still two leads:
(1) Maybe I'm doing something wrong with FME raster objects and GDAL objects,
(2) I'll try your suggestions and also DEMDistanceCalculator with a flat DEM raster?
+15- Safer
Hi @fiorinoh,
I wonder whether I understand the proximity functionality correctly. Check the attached image - would you like to get pixels representing distances from blue pixels to the closest red pixels? If so (the raster on the right and the inset), I have a pure FME solution for this, but it's a bit inefficient at this point, I am looking for ways to make it faster at the moment.
Dmitri
Hi @fiorinoh,
I wonder whether I understand the proximity functionality correctly. Check the attached image - would you like to get pixels representing distances from blue pixels to the closest red pixels? If so (the raster on the right and the inset), I have a pure FME solution for this, but it's a bit inefficient at this point, I am looking for ways to make it faster at the moment.
Dmitri
Hi @dmitribagh I would like to compute with FME transformers a proximity map as it is done with gdal_proximity. I attached below 2 images, one with red roads (the input GeoTIFF), and the other one showing the resulting proximity map. Thanks for your interest. Best regards, Humbert.
+15- Safer
Hi @fiorinoh,
I think, we can do it using standard FME tools. Here is how I did it - first, I generate a distance raster where each pixel represents the distance between the center pixel and all other pixels (I use a point cloud approach for this). Then, I turn all the road network pixels to points and pass my distance raster to them. Then, I mosaic all rasters together using minimal pixel values for overlapping rasters. As a last step, I deal with pixels beyond Distance parameter setting them to NoData. I tried to do this with the original distance raster, but somehow, it significantly affects performance (we'll investigate).
I used the road network of the city where I live. Once rasterized (1m resolution), it creates over 144,000 pixels representing roads.
My laptop generates the output in 1 minute if I set distance to 200 meters (400x400 distance raster), in 2.5 minutes with 400 meters (800x800 pixels), and under 4 minutes for 500 m distance (1000x1000 pixels). The final raster is about 8K pixels along each side.
I'll see whether we can compare it somehow with the GDAL proximity calculator. I also can make a custom transformer from this workspace if you think this method does what you expect and its performance is reasonable. See the attached template. It was made with FME 2022, but it opens and runs with FME 2021.2 as well.
Dmitri
+15- Safer
Hi @dmitribagh I would like to compute with FME transformers a proximity map as it is done with gdal_proximity. I attached below 2 images, one with red roads (the input GeoTIFF), and the other one showing the resulting proximity map. Thanks for your interest. Best regards, Humbert.
Hi @fiorinoh , see my post below. I am not sure you were properly tagged in it. Sorry about that.
Dmitri
+18- Supporter
You can kind of brute force it with a NeighbourFinder, but it might not be very efficient.
Get a bounding box around your data, then 2DGridAccumulator to convert it to a grid of points. NeighborFinder vs the road lines, 3DForcer based on _distance, and NumericRasterizer to create the final image.
If you want cell sizes based on your original raster, you can resample it down/clip it to a buffer of the roads, then RasterCellCoercer it into points to pass into the NeighborFinder.
Using Dmitri's road network sample with a 10m resolution on the 2DGridAccumulator, it's a little over a minute to process. Obviously going down to a 1m resolution will mean processing 100x the number of points, and so the mileage of this approach may vary.
Hi @fiorinoh,
I think, we can do it using standard FME tools. Here is how I did it - first, I generate a distance raster where each pixel represents the distance between the center pixel and all other pixels (I use a point cloud approach for this). Then, I turn all the road network pixels to points and pass my distance raster to them. Then, I mosaic all rasters together using minimal pixel values for overlapping rasters. As a last step, I deal with pixels beyond Distance parameter setting them to NoData. I tried to do this with the original distance raster, but somehow, it significantly affects performance (we'll investigate).
I used the road network of the city where I live. Once rasterized (1m resolution), it creates over 144,000 pixels representing roads.
My laptop generates the output in 1 minute if I set distance to 200 meters (400x400 distance raster), in 2.5 minutes with 400 meters (800x800 pixels), and under 4 minutes for 500 m distance (1000x1000 pixels). The final raster is about 8K pixels along each side.
I'll see whether we can compare it somehow with the GDAL proximity calculator. I also can make a custom transformer from this workspace if you think this method does what you expect and its performance is reasonable. See the attached template. It was made with FME 2022, but it opens and runs with FME 2021.2 as well.
Dmitri
Thanks a lot for this impressive workspace. Yes I think it would be great to have a transformer that does proximity calculation as GDAL. I'm also comparing your workspace with GDAL proximity calculator. BR.
You can kind of brute force it with a NeighbourFinder, but it might not be very efficient.
Get a bounding box around your data, then 2DGridAccumulator to convert it to a grid of points. NeighborFinder vs the road lines, 3DForcer based on _distance, and NumericRasterizer to create the final image.
If you want cell sizes based on your original raster, you can resample it down/clip it to a buffer of the roads, then RasterCellCoercer it into points to pass into the NeighborFinder.
Using Dmitri's road network sample with a 10m resolution on the 2DGridAccumulator, it's a little over a minute to process. Obviously going down to a 1m resolution will mean processing 100x the number of points, and so the mileage of this approach may vary.
I did some tests with @dmitribagh 's approach, and, as with your proposition, it seems the main issue is brute force and cache-memory. I'll try to figure out a workaround but, yes, it works!
Hi @fiorinoh,
I think, we can do it using standard FME tools. Here is how I did it - first, I generate a distance raster where each pixel represents the distance between the center pixel and all other pixels (I use a point cloud approach for this). Then, I turn all the road network pixels to points and pass my distance raster to them. Then, I mosaic all rasters together using minimal pixel values for overlapping rasters. As a last step, I deal with pixels beyond Distance parameter setting them to NoData. I tried to do this with the original distance raster, but somehow, it significantly affects performance (we'll investigate).
I used the road network of the city where I live. Once rasterized (1m resolution), it creates over 144,000 pixels representing roads.
My laptop generates the output in 1 minute if I set distance to 200 meters (400x400 distance raster), in 2.5 minutes with 400 meters (800x800 pixels), and under 4 minutes for 500 m distance (1000x1000 pixels). The final raster is about 8K pixels along each side.
I'll see whether we can compare it somehow with the GDAL proximity calculator. I also can make a custom transformer from this workspace if you think this method does what you expect and its performance is reasonable. See the attached template. It was made with FME 2022, but it opens and runs with FME 2021.2 as well.
Dmitri
Hi @dmitribagh , this is my last solution to this issue. As it is not easy to compute proximity rasters with FME transformers, I decided to use a python caller + GDAL. That was not easy because FME raster type and GDAL raster type are not compatible. Hope this code will help the FME community. BR.
import fmeobjects
import numpy as np
from osgeo import gdal
class MyBandTilePopulator(fmeobjects.FMEBandTilePopulator):
"""
This is a subclass of the FMEBandTilePopulator superclass.
It will be used when data is requested to create a new tile and
and populate it to a new FMEBand.
"""
def __init__(self, rasterData):
self._rasterData = rasterData
# required method
def clone(self):
"""
This method is used to create a copy of the data
multiple times while creating a new band
"""
return MyBandTilePopulator(self._rasterData)
# required method
def getTile(self, startRow, startCol, tile):
"""
Creates a new tile that's sized based on the input tile.
Populates that tile using this populator's raster data beginning
at the startRow and startCol.
"""
numRows, numCols = tile.getNumRows(), tile.getNumCols()
newTile = fmeobjects.FMEUInt16Tile(numRows, numCols)
data = newTile.getData()
for row in range(startRow, startRow+numRows):
for col in range(startCol, startCol+numCols):
if row < len(self._rasterData) and col < len(self._rasterData[0]):
data[row - startRow][col - startCol] = self._rasterData[row][col]
newTile.setData(data)
return newTile
class FeatureProcessor(object):
def __init__(self):
self.rasterData = []
def input(self, feature):
self.sysRef = feature.getCoordSys()
raster = feature.getGeometry()
rp = raster.getProperties()
band = raster.getBand(0)
bp = band.getProperties()
print("=== FME Input Raster ===")
print("Coordinate System = " + self.sysRef)
self.pixelWidth = rp.getSpacingX()
print("pixelWidth = " + str(self.pixelWidth))
self.pixelHeight = rp.getSpacingY()
print("pixelHeight = " + str(self.pixelHeight))
self.originX = rp.getOriginX()
print("originX = " + str(self.originX))
self.originY = rp.getOriginY()
print("originY = " + str(self.originY))
numRows = rp.getNumRows()
print("numRows = " + str(numRows))
numCols = rp.getNumCols()
print("numCols = " + str(numCols))
print("===")
tile = fmeobjects.FMEGray16Tile(numRows, numCols)
bandData = band.getTile(0, 0, tile).getData()
array = np.array(bandData)
driver = gdal.GetDriverByName('GTiff')
num_of_bands = 1
print("Creating file raster_transportation_roads.tiff")
srcRaster = driver.Create('raster_transportation_roads.tiff', numCols, numRows, num_of_bands, gdal.GDT_UInt16)
srcRaster.SetGeoTransform((self.originX, self.pixelWidth, 0, self.originY, 0, self.pixelHeight))
srcBand = srcRaster.GetRasterBand(1)
srcBand.WriteArray(array)
srcBand.FlushCache()
srcBand = None
srcRaster = None
src_ds = gdal.Open('raster_transportation_roads.tiff')
geotransform = src_ds.GetGeoTransform()
srcBand = src_ds.GetRasterBand(1)
print("Origin = ({}, {})".format(geotransform[0], geotransform[3]))
print("Pixel Size = ({}, {})".format(geotransform[1], geotransform[5]))
print("===")
print("Creating file proximity_transportation_roads.tiff with GDAL")
dst_filename='proximity_transportation_roads.tiff'
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create(dst_filename, numCols, numRows, num_of_bands, gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform(geotransform)
dst_ds.SetProjection(src_ds.GetProjectionRef())
dstBand = dst_ds.GetRasterBand(1)
gdal.ComputeProximity(srcBand, dstBand, ["DISTUNITS=PIXEL"])
self.rasterData = dstBand.ReadAsArray().tolist()
dstBand.FlushCache()
dstBand = None
srcBand = None
print("===")
def close(self):
# creating the raster data and specifying the formatting of the new raster
rasterData = self.rasterData
# specifying all of the properties for the new FMERaster
numRows, numCols = len(rasterData), len(rasterData[0])
xCellOrigin, yCellOrigin = 0.5, 0.5
xSpacing, ySpacing = self.pixelWidth, self.pixelHeight
xOrigin, yOrigin = self.originX, self.originY
xRotation, yRotation = 0.0, 0.0
# creating the new FMERaster
rasterProperties = fmeobjects.FMERasterProperties(numRows, numCols,
xSpacing, ySpacing,
xCellOrigin, yCellOrigin,
xOrigin, yOrigin,
xRotation, yRotation)
raster = fmeobjects.FMERaster(rasterProperties)
# Populating the contents of the band and appending it to the raster
bandTilePopulator = MyBandTilePopulator(rasterData)
bandName = ''
bandProperties = fmeobjects.FMEBandProperties(bandName,
fmeobjects.FME_INTERPRETATION_UINT16,
fmeobjects.FME_TILE_TYPE_FIXED,
numRows, numCols)
band = fmeobjects.FMEBand(bandTilePopulator, rasterProperties,
bandProperties)
raster.appendBand(band)
# creating a new feature with the FMERaster geometry to be output
feature = fmeobjects.FMEFeature()
feature.setGeometry(raster)
feature.setCoordSys(self.sysRef)
self.pyoutput(feature)
def process_group(self):
"""When 'Group By' attribute(s) are specified, this method is called
once all the FME Features in a current group have been sent to input().
FME Features sent to input() should generally be cached for group-by
processing in this method when knowledge of all Features is required.
The resulting Feature(s) from the group-by processing should be emitted
through self.pyoutput().
FME will continue calling input() a number of times followed
by process_group() for each 'Group By' attribute, so this
implementation should reset any class members for the next group.
"""
pass
Reply
Helpful Members This Week
Enter your E-mail address. We'll send you an e-mail with instructions to reset your password.