JuliaDataCubes / PyramidScheme.jl

Building and using pyramids for large raster data
MIT License
10 stars 1 forks source link

Add plot dispatch on GeoAxis and Pyramid #45

Open felixcremer opened 1 week ago

felixcremer commented 1 week ago

If we would like to combine the pyramid data with data in other projections it would be best when the pyramid data is reprojected on the fly during interactions to the target projection of the plotting Axis. Therefore we would like to add this interaction into a dispatch of GeoAxis. The code for the reprojection on the fly is

    on(ax.finallimits) do limits
        limext = Extents.extent(limits)
        # Compute limit in raster projection
        #todataproj = Proj.Transformation(plotcrs, rastercrs)
        #@show plotcrs, rastercrs
        trans = Proj.Transformation(plotcrs, rastercrs, always_xy=true)
        datalimit = trans_bounds(trans, limext)
        #datalimit = limext
        #@show datalimit
        if Extents.intersects(rasext, datalimit)
            rasdata = selectlevel(pyramids, datalimit)
            # Project selected data to plotcrs
            #data.val = rasdata
            data.val = Rasters.resample(rasdata, crs="EPSG:3857")
        end
        notify(data)
    end

This is from this commit: https://github.com/JuliaDataCubes/PyramidScheme.jl/blob/fb8eba07376f04aed3b424e2994321fee0b58f0f/src/PyramidScheme.jl

Parts of this are still in the source code and should be taken out.

asinghvi17 commented 1 week ago

The code I see above won't work on general transformed axes (GeoAxis included) since you will pick up the bounding box in transformed space. Transforming back to input space is also a lossy operation - consider an axis aligned square centered at 0, now rotate it by 45 degrees. The bounding box of this square has double the area of the original square.

This shouldn't be an issue for interactive use but you have to be careful about how you specify the space and transformations.

felixcremer commented 1 week ago

Also I am currently use the viewport widths as indicator how large the data should be that im grabbing. This might be wrong for geoaxies.