locationtech / rasterframes

Geospatial Raster support for Spark DataFrames
http://rasterframes.io
Apache License 2.0
240 stars 46 forks source link

Request new function rf_agg_local_percentile #575

Open lepitrust opened 2 years ago

lepitrust commented 2 years ago

The function proto should be:

Tile rf_agg_local_percentile(Tile tile, percentile, List[float] probabilities, float relative_error)

The result shoud be a tile with same crs - extent of input tiles (same way of rf_agg_local_max etc), eventually the probability could be only one and not a list. Percentile shoud be calculate cell-wise over tile column.

https://docs.qgis.org/3.22/en/docs/user_manual/processing_algs/qgis/rasteranalysis.html#cell-stack-percentile

Thanks

EDIT: In addiction to anyone do the same thing now i write a code, the column "ghi" is a tile column. The ghi are multiple raster chunked by rasterframes to tile 256x256. Group by extent and crs to maintain align over the tile. Use xmin ymin etc to speed up the goupby. The dims is necessary to rebuild correctly the tile on output.

df_rf_ghi = df_rf_ghi.withColumn('extent', st_extent(rf_geometry('ghi')))
df_rf_ghi = df_rf_ghi.withColumn('dims', rf_dimensions('ghi'))
df_rf_ghi = df_rf_ghi.withColumn('crs', rf_crs('ghi'))
df_rf_ghi = df_rf_ghi.select(F.col('my_index'), F.col('extent'), F.col('crs'), F.col('extent')["xmin"].alias("xmin"), \
                 F.col('extent')["ymin"].alias("ymin"), F.col('extent')["xmax"].alias("xmax"), F.col('extent')["ymax"].alias("ymax"), \
                 F.col('dims')["cols"].alias("cols"), F.col('dims')["rows"].alias("rows"), rf_explode_tiles('ghi'))
med = df_rf_ghi.dropna(subset="ghi")
med = med.groupby('xmin', 'xmax', 'ymin', 'ymax', 'cols', 'rows', 'extent', 'crs', 'column_index', 'row_index') \
            .agg(F.expr('percentile_approx(ghi, array(0.75), 100)')[0].alias('_75'))
med.createOrReplaceTempView("med");
rebuild = spark.sql("SELECT extent, crs, rf_assemble_tile(column_index, row_index, _75 , cols, rows) AS tile_percentile \
            FROM med GROUP BY xmin, xmax, ymin, ymax, cols, rows, extent, crs")
rebuild2 = rebuild.select(rf_proj_raster(F.col('tile_percentile'), F.col('extent'), F.col('crs')).alias('pro_ghi'))