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?

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.

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 @fiorinoh​ , see my post below. I am not sure you were properly tagged in it. Sorry about that.

 

Dmitri

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.

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.

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 @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._rasterDatai0]):
                   datatrow - 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(geotransformD0], geotransformT3]))
        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, s"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