[BUG] gridding_plot only work for global and WGS84 currently

[juangallegozamorano]

Hi @chenyangkang,
Exporting the plot for the grid of the quadtrees only works when using WGS84. Is there a way to access the grid itself and or export the figure under another CRS?

Thanks!

[chenyangkang]

Hi @juangallegozamorano !
Thanks for opening the issue!

1. The grid plot extent is currently set to (-180,180) and (-90, 90). You can changed it to your CRS extent by setting plot_xlims and plot_ylims params during model class initiation. For example plot_xlims=(-64000, 64000). I will add an auto-estimator in the next PR.
2. The girds are generated based on your input data. So it makes less sense to project the grids to another CRS, and imaginably they will be twisted if you do that.
3. You can still access the grids by calling model.ensemble_df.

[juangallegozamorano]

Hi @chenyangkang, I accessed the ensemble_df but unfortunately I'm not able to plot the stixels. Can you tell me how to do it? Thanks!

[chenyangkang]

@juangallegozamorano There is not function implemented for this in stemflow. But you can do:

from stemflow.utils.jitterrotation.jitterrotator import JitterRotator
from shapely.geometry import Polygon
import geopandas as gpd

# define a function
def geo_grid_geometry(line):
    old_x, old_y = JitterRotator.inverse_jitter_rotate(
        [line['stixel_calibration_point_transformed_left_bound'], line['stixel_calibration_point_transformed_left_bound'], line['stixel_calibration_point_transformed_right_bound'], line['stixel_calibration_point_transformed_right_bound']],
        [line['stixel_calibration_point_transformed_lower_bound'], line['stixel_calibration_point_transformed_upper_bound'], line['stixel_calibration_point_transformed_upper_bound'], line['stixel_calibration_point_transformed_lower_bound']],
        line['rotation'],
        line['calibration_point_x_jitter'],
        line['calibration_point_y_jitter'],
    )

    polygon = Polygon(list(zip(old_x, old_y)))
    return polygon

# Make a geometry attribute for each stixel
model.ensemble_df['geometry'] = model.ensemble_df.apply(geo_grid_geometry, axis=1)
model.ensemble_df = gpd.GeoDataFrame(model.ensemble_df, geometry='geometry')

Which creates a spatial object for each stixel.

Then:

model.ensemble_df.plot(alpha=0.2)

The stixel will stack on the temporal dimension. You may want to pick only one time range (since the spatial splitting pattern will be the same for all temporal windows in one ensemble).

model.ensemble_df[
        (model.ensemble_df['DOY_start']>=90) & (model.ensemble_df['DOY_start']<120)
].plot(alpha=0.2)

Which will look less messy.

You may also want to add scatter points on top of this:

fig,ax = plt.subplots()

model.ensemble_df[
        (model.ensemble_df['DOY_start']>=90) & ((model.ensemble_df['DOY_start']<120))
].plot(alpha=0.2,ax=ax)

ax.scatter(
        X_train['longitude'],
        X_train['latitude'],
        c='tab:orange',s=0.2,alpha=0.7
)

plt.show()

Cheers!

[juangallegozamorano]

Thanks a lot!! This is really helpful!

[juangallegozamorano]

Hi @chenyangkang, for some reason the geometries are all empty in my case...

[chenyangkang]

@juangallegozamorano You could provide details e.g. the model.ensemble_df original data frame for debugging purpose.

[juangallegozamorano]

Here are two examples of ensemble_df that I would like to represent in a figure
53970_larcan_ensemble_grid.csv
Common_Gull_ensemble_df.csv

[chenyangkang]

@juangallegozamorano I cannot reproduce the bug. This is what I got for 53970_larcan_ensemble_grid.csv:

data = pd.read_csv('53970_larcan_ensemble_grid.csv', index_col=0)

from stemflow.utils.jitterrotation.jitterrotator import JitterRotator
from shapely.geometry import Polygon
import geopandas as gpd

# define a function
def geo_grid_geometry(line):
    old_x, old_y = JitterRotator.inverse_jitter_rotate(
        [line['stixel_calibration_point_transformed_left_bound'], line['stixel_calibration_point_transformed_left_bound'], line['stixel_calibration_point_transformed_right_bound'], line['stixel_calibration_point_transformed_right_bound']],
        [line['stixel_calibration_point_transformed_lower_bound'], line['stixel_calibration_point_transformed_upper_bound'], line['stixel_calibration_point_transformed_upper_bound'], line['stixel_calibration_point_transformed_lower_bound']],
        line['rotation'],
        line['calibration_point_x_jitter'],
        line['calibration_point_y_jitter'],
    )

    polygon = Polygon(list(zip(old_x, old_y)))
    return polygon

# Make a geometry attribute for each stixel
data['geometry'] = data.apply(geo_grid_geometry, axis=1)
data = gpd.GeoDataFrame(data, geometry='geometry')

[chenyangkang]

I also successfully plotted the Common_Gull_ensemble_df.csv. Maybe you could provide more information? For example:

add a "print" line in the function:

def geo_grid_geometry(line):
    old_x, old_y = JitterRotator.inverse_jitter_rotate(
        [line['stixel_calibration_point_transformed_left_bound'], line['stixel_calibration_point_transformed_left_bound'], line['stixel_calibration_point_transformed_right_bound'], line['stixel_calibration_point_transformed_right_bound']],
        [line['stixel_calibration_point_transformed_lower_bound'], line['stixel_calibration_point_transformed_upper_bound'], line['stixel_calibration_point_transformed_upper_bound'], line['stixel_calibration_point_transformed_lower_bound']],
        line['rotation'],
        line['calibration_point_x_jitter'],
        line['calibration_point_y_jitter'],
    )

    polygon = Polygon(list(zip(old_x, old_y)))
    print(polygon) # this line should print out the polygon string.
    return polygon

[juangallegozamorano]

I just tried again and now I can see the stixels...I cannot explain why there were empty geometries before, thanks a lot and sorry for the hassle!

[chenyangkang]

No worries! Good to hear you solved it!