Question

Calculate the median of overlapping pixel values from 24 input raster bands

U
Badge

Dear FME users,

 

 

I’d like to calculate the median of all overlapping pixel values from 24 input raster bands.

 

 

Do you have an idea how to solve it? Can I use the transformer “RasterExpressionEvaluator”?

 

 

Thank you very much.

 

 

Konrad

21 replies

G
Badge +3

expressionevaluator has only A en B input.

B can have multiple raster but will be paired, so B's do not calculate between them.

You could use a wsc to sequentially calculate all B inputs with a updating A input.

Or first extract the overlapping bits, extract the data, process the data and rebuild raster.

jdh
Rank
Badge +22

The last time I did something like that was several years ago, where I coerced the rasters to points, aggregated the points on col, row keeping a list, and then sorted and analysed the list.

Note that this was for 24 1 band rasters, if you have a single raster with 24 bands, you can coerce the rasterer with bands as attributes and skip the aggregation step.

 

If I were to do it now, I would first look into coercing to pointclouds, and if that didn't lead anywhere, then experiment with the new python raster api.

Note that if you wanted the mean rather than the medium, then the RasterMosaicker with Overlapping values set to Average would be the way to go.

U
Badge

Thank you for all your inputs @jdh and @gio.

 

 

I’ve tried to calculate it with the transformer “RasterExpressionEvaluator”. I think it works, but the expressions are overloaded.

 

 

Here is the first line of the repeated code:

 

“@if(A[23],(A[0]+A[1]+A[2]+A[3]+A[4]+A[5]+A[6]+A[7]+A[8]+A[9]+A[10]+A[11]+A[12]+A[13]+A[14]+A[15]+A[16]+A[17]+A[18]+A[19]+A[20]+A[21]+A[22]+A[23])/24,A[0])”

 

 

Are there any ideas how to simplify the expressions? It's possible to have just one expression?

 

 

I’ve added some screenshots in the attachment.

 

 

Thank you very much. :)

 

 

Best regards. Konrad
david_r
Rank
Userlevel 4

The RasterCellCoercer mentioned by @jdh can help you along, just be aware that it can be very slow when dealing with large rasters. Apparently you can now coerce the raster to a point cloud, which should be a lot faster, but I've got no experience with that.

If you're no stranger to Python, I'd recommend looking into the new raster capabilities of the fmeobjects API. @takashi wrote some really nice sample code here, which could be a starting point:

https://knowledge.safe.com/questions/38000/python-fme-objects-api-for-raster-manipulation.html

david_r
Rank
Userlevel 4

Thank you for all your inputs @jdh and @gio.

 

 

I’ve tried to calculate it with the transformer “RasterExpressionEvaluator”. I think it works, but the expressions are overloaded.

 

 

Here is the first line of the repeated code:

 

“@if(A[23],(A[0]+A[1]+A[2]+A[3]+A[4]+A[5]+A[6]+A[7]+A[8]+A[9]+A[10]+A[11]+A[12]+A[13]+A[14]+A[15]+A[16]+A[17]+A[18]+A[19]+A[20]+A[21]+A[22]+A[23])/24,A[0])”

 

 

Are there any ideas how to simplify the expressions? It's possible to have just one expression?

 

 

I’ve added some screenshots in the attachment.

 

 

Thank you very much. :)

 

 

Best regards. Konrad
I'm not sure if it's possible to simplify it much more. The RasterExpressionEvaluator is a bit old-fashioned and quirky in the way it works, but if it suits your use-case it is pretty fast.
takashi
Rank
Userlevel 2
Badge +17

Hi @UBA_KP, if you need to calculate the "median" of 24 band values for each cell, the RasterExpressionEvaluator won't help you. Because, the transformer doesn't support list (array) manipulation, although it's necessary to create a list containing the 24 values and sort the list in order to compute the "median". e.g.

List: [1, 3, 1, 3, 2, 1] -> Sorted: [1, 1, 1, 2, 3, 3] -> Median: (1 + 2) / 2 = 1.5

See here to lean more about "median": https://en.wikipedia.org/wiki/Median

As @jdh and @david_r suggested, you can get individual grid points each of them has 24 band values as a list attribute with coercing the raster with the RasterCellCoercer. You can then use the ListSorter and the ExpressionEvaluator to calculate the median base on the list. Finally set the median value to z-coordinate to each point with the 3DForcer and re-create a raster with a NumericRasterizer.

However, the performance of coercing is not good. If you don't like the coercing approach, Python scripting could be a better solution, as @david_r mentioned.

Alternatively, this could also be a possible way.

  1. RasterBandSeparator: Separate the source raster into 24 rasters each of them has a single band.
  2. PointCloudCombiner: Create a single point cloud from the 24 rasters. The point cloud will have 24 points for every grid cell.
  3. PointCloudStatsRasterizer (from Hub): This transformer calculates a statistics value based on the multiple points for every grid cell covering the input point cloud. Set "median" to the Statistics Calculation Mode parameter and set the X/Y cell spacing of the source raster to the X/Y Cell Spacing parameters. You can use the RasterPropertyExtractor to extract the X/Y cell spacing of the source raster.
U
Badge

Hi @UBA_KP, if you need to calculate the "median" of 24 band values for each cell, the RasterExpressionEvaluator won't help you. Because, the transformer doesn't support list (array) manipulation, although it's necessary to create a list containing the 24 values and sort the list in order to compute the "median". e.g.

List: [1, 3, 1, 3, 2, 1] -> Sorted: [1, 1, 1, 2, 3, 3] -> Median: (1 + 2) / 2 = 1.5

See here to lean more about "median": https://en.wikipedia.org/wiki/Median

As @jdh and @david_r suggested, you can get individual grid points each of them has 24 band values as a list attribute with coercing the raster with the RasterCellCoercer. You can then use the ListSorter and the ExpressionEvaluator to calculate the median base on the list. Finally set the median value to z-coordinate to each point with the 3DForcer and re-create a raster with a NumericRasterizer.

However, the performance of coercing is not good. If you don't like the coercing approach, Python scripting could be a better solution, as @david_r mentioned.

Alternatively, this could also be a possible way.

  1. RasterBandSeparator: Separate the source raster into 24 rasters each of them has a single band.
  2. PointCloudCombiner: Create a single point cloud from the 24 rasters. The point cloud will have 24 points for every grid cell.
  3. PointCloudStatsRasterizer (from Hub): This transformer calculates a statistics value based on the multiple points for every grid cell covering the input point cloud. Set "median" to the Statistics Calculation Mode parameter and set the X/Y cell spacing of the source raster to the X/Y Cell Spacing parameters. You can use the RasterPropertyExtractor to extract the X/Y cell spacing of the source raster.
Thank your very much @takashi. Maybe you can give me an example of your first recommendation (list). I gues that I'll need following transformers: ListBuilder, ListSorter and StatisticsCalculator. I've forgotten to tell you that I use the transformer RasterToPolygonCoercer before.

 

takashi
Rank
Userlevel 2
Badge +17

Hi @UBA_KP, if you need to calculate the "median" of 24 band values for each cell, the RasterExpressionEvaluator won't help you. Because, the transformer doesn't support list (array) manipulation, although it's necessary to create a list containing the 24 values and sort the list in order to compute the "median". e.g.

List: [1, 3, 1, 3, 2, 1] -> Sorted: [1, 1, 1, 2, 3, 3] -> Median: (1 + 2) / 2 = 1.5

See here to lean more about "median": https://en.wikipedia.org/wiki/Median

As @jdh and @david_r suggested, you can get individual grid points each of them has 24 band values as a list attribute with coercing the raster with the RasterCellCoercer. You can then use the ListSorter and the ExpressionEvaluator to calculate the median base on the list. Finally set the median value to z-coordinate to each point with the 3DForcer and re-create a raster with a NumericRasterizer.

However, the performance of coercing is not good. If you don't like the coercing approach, Python scripting could be a better solution, as @david_r mentioned.

Alternatively, this could also be a possible way.

  1. RasterBandSeparator: Separate the source raster into 24 rasters each of them has a single band.
  2. PointCloudCombiner: Create a single point cloud from the 24 rasters. The point cloud will have 24 points for every grid cell.
  3. PointCloudStatsRasterizer (from Hub): This transformer calculates a statistics value based on the multiple points for every grid cell covering the input point cloud. Set "median" to the Statistics Calculation Mode parameter and set the X/Y cell spacing of the source raster to the X/Y Cell Spacing parameters. You can use the RasterPropertyExtractor to extract the X/Y cell spacing of the source raster.
Since the raster has 24 bands, the RasterCellCoercer (Extract Band Values: Attributes) adds a list attribute called "_band{}.value" to every resulting point, and the list contains 24 elements storing the band values of the source cell. You can just sort the list and then calculate the median. This screenshot illustrates the procedure.

 

 

U
Badge
Since the raster has 24 bands, the RasterCellCoercer (Extract Band Values: Attributes) adds a list attribute called "_band{}.value" to every resulting point, and the list contains 24 elements storing the band values of the source cell. You can just sort the list and then calculate the median. This screenshot illustrates the procedure.

 

 

Thank you so much @takashi. I've done a mistake in my question. I'll desribe it very soon.

 

U
Badge

Dear FME-Users,

 

 

I’ve done a mistake during the translation. Sorry, I’m german. :))) What I need is not the median but the running/moving average is what I’d like to calculate.

 

 

Now I repeat myself again:

 

I’ve tried to calculate it with the transformer RasterExpressionEvaluator. I think it works, but the expressions are overloaded because all raster have 73 bands.

 

 

Here is the first line from my code:

 

“@if(A[23],(A[0]+A[1]+A[2]+A[3]+A[4]+A[5]+A[6]+A[7]+A[8]+A[9]+A[10]+A[11]+A[12]+A[13]+A[14]+A[15]+A[16]+A[17]+A[18]+A[19]+A[20]+A[21]+A[22]+A[23])/24,A[0])”

 

Here is the last line from my code:

 

“@if(A[73],(A[50]+A[51]+A[52]+A[53]+A[54]+A[55]+A[56]+A[57]+A[58]+A[59]+A[60]+A[61]+A[62]+A[63]+A[64]+A[65]+A[66]+A[67]+A[68]+A[69]+A[70]+A[71]+A[72]+A[73])/24,A[50])”

 

 

Are there any ideas how to simplify the expressions? It's possible to have just one expression?

 

 

After the transformer RasterExpressionEvaluator I’ve connected the transformer RasterToPolygonCoercer to generate vector data.

 

 

Thank you very much and best regards!

 

 

Konrad
takashi
Rank
Userlevel 2
Badge +17

moving average...

If you are using FME 2017.0+, the RasterMosaicker can be used to calculate the average of multiple band values, as @jdh suggested before. A possible way to get moving average in 73 bands for each moving range consisting of 24 bands is:

  1. Cloner (Number of Copies: 50, Copy Number Attribute: _copynum): Create 50 (= 73 - 24 + 1) copies of the source raster. 50 is the number of moving ranges.
  2. RasterBandSeparator (Band Index Attribute: _band_index): Decompose the copied rasters into single band rasters. [50 x 73 = 3,650 single band rasters]
  3. Tester (_copy_num <= _band_index AND _band_index < _copy_num + 24): Select rasters for each moving range (24 bands). [50 x 24 = 1,200 single band rasters]
  4. RasterMosaicker (Group By: _copynum, Input Ordered: By Group, Overlapping Value: Average): Mosaick rasters for each moving range (24 bands). Here, the average of 24 bands will be resulting cell values if you set Average to the Overlapping Value parameter. [50 single band rasters]
  5. RasterBandCombiner (Group By: <not set>): Combine the rasters into a single raster having 50 bands, which contain moving average values. [a single raster consisting of 50 bands]

However, I cannot understand why the @if function appears in your expression examples. I don't think conditional operation is necessary to calculate average. If there is a specific reason for that a conditional operation is required, the procedure above should be modified (or could not be applied).

U
Badge

moving average...

If you are using FME 2017.0+, the RasterMosaicker can be used to calculate the average of multiple band values, as @jdh suggested before. A possible way to get moving average in 73 bands for each moving range consisting of 24 bands is:

  1. Cloner (Number of Copies: 50, Copy Number Attribute: _copynum): Create 50 (= 73 - 24 + 1) copies of the source raster. 50 is the number of moving ranges.
  2. RasterBandSeparator (Band Index Attribute: _band_index): Decompose the copied rasters into single band rasters. [50 x 73 = 3,650 single band rasters]
  3. Tester (_copy_num <= _band_index AND _band_index < _copy_num + 24): Select rasters for each moving range (24 bands). [50 x 24 = 1,200 single band rasters]
  4. RasterMosaicker (Group By: _copynum, Input Ordered: By Group, Overlapping Value: Average): Mosaick rasters for each moving range (24 bands). Here, the average of 24 bands will be resulting cell values if you set Average to the Overlapping Value parameter. [50 single band rasters]
  5. RasterBandCombiner (Group By: <not set>): Combine the rasters into a single raster having 50 bands, which contain moving average values. [a single raster consisting of 50 bands]

However, I cannot understand why the @if function appears in your expression examples. I don't think conditional operation is necessary to calculate average. If there is a specific reason for that a conditional operation is required, the procedure above should be modified (or could not be applied).

Thank you very much @takashi! But I don't have FME 2017.0+, our version is FME 2016.1. Which transformer can I use instead of RasterMosaicker?

 

takashi
Rank
Userlevel 2
Badge +17

moving average...

If you are using FME 2017.0+, the RasterMosaicker can be used to calculate the average of multiple band values, as @jdh suggested before. A possible way to get moving average in 73 bands for each moving range consisting of 24 bands is:

  1. Cloner (Number of Copies: 50, Copy Number Attribute: _copynum): Create 50 (= 73 - 24 + 1) copies of the source raster. 50 is the number of moving ranges.
  2. RasterBandSeparator (Band Index Attribute: _band_index): Decompose the copied rasters into single band rasters. [50 x 73 = 3,650 single band rasters]
  3. Tester (_copy_num <= _band_index AND _band_index < _copy_num + 24): Select rasters for each moving range (24 bands). [50 x 24 = 1,200 single band rasters]
  4. RasterMosaicker (Group By: _copynum, Input Ordered: By Group, Overlapping Value: Average): Mosaick rasters for each moving range (24 bands). Here, the average of 24 bands will be resulting cell values if you set Average to the Overlapping Value parameter. [50 single band rasters]
  5. RasterBandCombiner (Group By: <not set>): Combine the rasters into a single raster having 50 bands, which contain moving average values. [a single raster consisting of 50 bands]

However, I cannot understand why the @if function appears in your expression examples. I don't think conditional operation is necessary to calculate average. If there is a specific reason for that a conditional operation is required, the procedure above should be modified (or could not be applied).

Correction.

 

  • RasterMosaicker: Set the Input Ordered parameter.
  • RasterBandCombiner: Not set the Group By parameter

 

takashi
Rank
Userlevel 2
Badge +17

moving average...

If you are using FME 2017.0+, the RasterMosaicker can be used to calculate the average of multiple band values, as @jdh suggested before. A possible way to get moving average in 73 bands for each moving range consisting of 24 bands is:

  1. Cloner (Number of Copies: 50, Copy Number Attribute: _copynum): Create 50 (= 73 - 24 + 1) copies of the source raster. 50 is the number of moving ranges.
  2. RasterBandSeparator (Band Index Attribute: _band_index): Decompose the copied rasters into single band rasters. [50 x 73 = 3,650 single band rasters]
  3. Tester (_copy_num <= _band_index AND _band_index < _copy_num + 24): Select rasters for each moving range (24 bands). [50 x 24 = 1,200 single band rasters]
  4. RasterMosaicker (Group By: _copynum, Input Ordered: By Group, Overlapping Value: Average): Mosaick rasters for each moving range (24 bands). Here, the average of 24 bands will be resulting cell values if you set Average to the Overlapping Value parameter. [50 single band rasters]
  5. RasterBandCombiner (Group By: <not set>): Combine the rasters into a single raster having 50 bands, which contain moving average values. [a single raster consisting of 50 bands]

However, I cannot understand why the @if function appears in your expression examples. I don't think conditional operation is necessary to calculate average. If there is a specific reason for that a conditional operation is required, the procedure above should be modified (or could not be applied).

The workflow I suggested above transform a raster containing 73 bands into a raster containing 50 bands, and each band contains the moving average of original 24 bands. i.e.

 

band 0 contains the average of original 0 to 23 band,

 

band 1 contains the average of original 1 to 24 band,

 

band 2 contains the average of original 2 to 25 band,

 

...

 

and band 49 contains the average of original 49 to 72 band.

 

Is it your desired result theoretically?

 

 

takashi
Rank
Userlevel 2
Badge +17

moving average...

If you are using FME 2017.0+, the RasterMosaicker can be used to calculate the average of multiple band values, as @jdh suggested before. A possible way to get moving average in 73 bands for each moving range consisting of 24 bands is:

  1. Cloner (Number of Copies: 50, Copy Number Attribute: _copynum): Create 50 (= 73 - 24 + 1) copies of the source raster. 50 is the number of moving ranges.
  2. RasterBandSeparator (Band Index Attribute: _band_index): Decompose the copied rasters into single band rasters. [50 x 73 = 3,650 single band rasters]
  3. Tester (_copy_num <= _band_index AND _band_index < _copy_num + 24): Select rasters for each moving range (24 bands). [50 x 24 = 1,200 single band rasters]
  4. RasterMosaicker (Group By: _copynum, Input Ordered: By Group, Overlapping Value: Average): Mosaick rasters for each moving range (24 bands). Here, the average of 24 bands will be resulting cell values if you set Average to the Overlapping Value parameter. [50 single band rasters]
  5. RasterBandCombiner (Group By: <not set>): Combine the rasters into a single raster having 50 bands, which contain moving average values. [a single raster consisting of 50 bands]

However, I cannot understand why the @if function appears in your expression examples. I don't think conditional operation is necessary to calculate average. If there is a specific reason for that a conditional operation is required, the procedure above should be modified (or could not be applied).

I found the RasterMosaicker (Step 4) in the procedure above can be replaced with RasterCombiner + RasterExpressionEvaluator.

 

RasterCombiner

 

  • Group By: _copynum
  • Input Ordered: By Group
RasterExpressionEvaluator (transform each raster consisting of 24 bands into single band raster storing the average of the 24 band values)

 

  • Expression: (A[0]+A[1]+A[2]+A[3]+A[4]+A[5]+A[6]+A[7]+A[8]+A[9]+A[10]+A[11]+A[12]+A[13]+A[14]+A[15]+A[16]+A[17]+A[18]+A[19]+A[20]+A[21]+A[22]+A[23])/24.0
U
Badge
I found the RasterMosaicker (Step 4) in the procedure above can be replaced with RasterCombiner + RasterExpressionEvaluator.

 

RasterCombiner

 

  • Group By: _copynum
  • Input Ordered: By Group
RasterExpressionEvaluator (transform each raster consisting of 24 bands into single band raster storing the average of the 24 band values)

 

  • Expression: (A[0]+A[1]+A[2]+A[3]+A[4]+A[5]+A[6]+A[7]+A[8]+A[9]+A[10]+A[11]+A[12]+A[13]+A[14]+A[15]+A[16]+A[17]+A[18]+A[19]+A[20]+A[21]+A[22]+A[23])/24.0
Yes, I've tried a similar solution. It works, thank you very much again @takashi!

 

U
Badge
The workflow I suggested above transform a raster containing 73 bands into a raster containing 50 bands, and each band contains the moving average of original 24 bands. i.e.

 

band 0 contains the average of original 0 to 23 band,

 

band 1 contains the average of original 1 to 24 band,

 

band 2 contains the average of original 2 to 25 band,

 

...

 

and band 49 contains the average of original 49 to 72 band.

 

Is it your desired result theoretically?

 

 

At the beginning I've worked with 24 bands of raster. Now there're 97 bands.

 

 

You're right, each band contains the moving average of other 24 bands. It means that I've to transform a raster containing 97 bands into a raster with includes 73 bands.

 

 

Thank you very much @takashi.
takashi
Rank
Userlevel 2
Badge +17

moving average...

If you are using FME 2017.0+, the RasterMosaicker can be used to calculate the average of multiple band values, as @jdh suggested before. A possible way to get moving average in 73 bands for each moving range consisting of 24 bands is:

  1. Cloner (Number of Copies: 50, Copy Number Attribute: _copynum): Create 50 (= 73 - 24 + 1) copies of the source raster. 50 is the number of moving ranges.
  2. RasterBandSeparator (Band Index Attribute: _band_index): Decompose the copied rasters into single band rasters. [50 x 73 = 3,650 single band rasters]
  3. Tester (_copy_num <= _band_index AND _band_index < _copy_num + 24): Select rasters for each moving range (24 bands). [50 x 24 = 1,200 single band rasters]
  4. RasterMosaicker (Group By: _copynum, Input Ordered: By Group, Overlapping Value: Average): Mosaick rasters for each moving range (24 bands). Here, the average of 24 bands will be resulting cell values if you set Average to the Overlapping Value parameter. [50 single band rasters]
  5. RasterBandCombiner (Group By: <not set>): Combine the rasters into a single raster having 50 bands, which contain moving average values. [a single raster consisting of 50 bands]

However, I cannot understand why the @if function appears in your expression examples. I don't think conditional operation is necessary to calculate average. If there is a specific reason for that a conditional operation is required, the procedure above should be modified (or could not be applied).

If the source dataset consists of multiple single-band rasters (not a single raster containing multiple bands), the workflow I suggested above should be modified like this (FME 2016.1).

 

 

 

 

U
Badge
If the source dataset consists of multiple single-band rasters (not a single raster containing multiple bands), the workflow I suggested above should be modified like this (FME 2016.1).

 

 

 

 

Thank you very much @takashi! Sorry for my late answer, I was on holiday and afterwards I’ve done other things.

 

 

I’ve added and connected the transformers “RasterBandCombiner”, “RasterExpressionEvaluator” and “RasterBandSeparator” because I’d like to calculate the moving average. Usually I've 73 input raster, for calculating the moving average it’s necessary to use 96 input raster. I’ve added 73 expressions for the “RasterExpressionEvaluator”, it’s a lot but it works. When I use the “RasterBandCombiner” I loose important attribute values. I know it’s possible to write attributes into a list, but I don’t know how to allocate them again. I also don’t know how to extract the last 73 attribute values from 96 values altogether. The background is that I need all names of 73 input data as string. It works fine without the transformers for the moving average. I’ve added a screenshot in the attachment for better comprehension.

 

 

Is there anyone who can help me?

 

 

Thank you and best regards!

 

 

Konrad

 

takashi
Rank
Userlevel 2
Badge +17
If the source dataset consists of multiple single-band rasters (not a single raster containing multiple bands), the workflow I suggested above should be modified like this (FME 2016.1).

 

 

 

 

A possible way is:

 

  1. Branch the data flow output from the RasterSubsetter into two streams. Connect one stream to the current workflow, i.e. the RasterBandCombiner.
  2. On the second stream, put a Counter to add a sequential number starting with -23.
  3. Add a FeatureMerter to the workspace. Send the features output from the RasterBandSeparator to the Requestor port, send the features output from the Counter on the second stream to the Supprier port, and then merge the supplier to the requestor on Requestor: "_band_index" (0-based sequential number) and Supplier: "_count" (sequential number starting with -23).
U
Badge
A possible way is:

 

  1. Branch the data flow output from the RasterSubsetter into two streams. Connect one stream to the current workflow, i.e. the RasterBandCombiner.
  2. On the second stream, put a Counter to add a sequential number starting with -23.
  3. Add a FeatureMerter to the workspace. Send the features output from the RasterBandSeparator to the Requestor port, send the features output from the Counter on the second stream to the Supprier port, and then merge the supplier to the requestor on Requestor: "_band_index" (0-based sequential number) and Supplier: "_count" (sequential number starting with -23).
Thank you very much @takashi. Now it works! I saw that you already wrote an answer to my question there: How to calculate the time for NetCDF-files?. I am sorry for this.

 

 

Best regards!

 

 

Konrad

 

Reply


Community Stats

28,670
Posts
109,365
Replies
36,531
Members
Terms & ConditionsCookie settings
A product of
Safe Software respectfully acknowledges that we live, learn and work on the traditional and unceded territories of the Kwantlen, Katzie, and Semiahmoo First Nations.