[jira] [Commented] (SIS-514) GridGeometry.subgrid(…) should accept an area of interest with less dimensions than the base grid

["Johann Sorel (Jira)" <ji...@apache.org>]

    [ https://issues.apache.org/jira/browse/SIS-514?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17369335#comment-17369335 ] 

Johann Sorel commented on SIS-514:
----------------------------------

Here are some functions I made which can bypasse the problem by spliting the intersection on each single crs of the resource GridGeometry 


{code:java}
    private GridGeometry intersect(GridGeometry resource, GridGeometry query, GridRoundingMode mode) throws FactoryException {
        try {
            return resource.derive().subgrid(query).build();
        } catch (Exception ex) {
            //do nothing
        }

        //fallback solution
        final CoordinateReferenceSystem resourceCrs = resource.getCoordinateReferenceSystem();
        final CoordinateReferenceSystem queryCrs = query.getCoordinateReferenceSystem();

        //find compatible axis
        final List<SingleCRS> resourceCrsComps = CRS.getSingleComponents(resourceCrs);
        final List<SingleCRS> queryCrsComps = CRS.getSingleComponents(queryCrs);
        final int[] resourceCrsCompsOffset = new int[resourceCrsComps.size()];
        final int[] queryCrsCompsOffset = new int[resourceCrsComps.size()];
        for (int i = 0, k = 0, n = resourceCrsComps.size(); i < n; i++) {
            final SingleCRS s = resourceCrsComps.get(i);
            resourceCrsCompsOffset[i] = k;
            k+= s.getCoordinateSystem().getDimension();
        }
        for (int i = 0, k = 0, n = queryCrsComps.size(); i < n; i++) {
            final SingleCRS s = queryCrsComps.get(i);
            queryCrsCompsOffset[i] = k;
            k+= s.getCoordinateSystem().getDimension();
        }

        final Map<SingleCRS,SingleCRS> map = new HashMap<>();
        final Map<SingleCRS,SingleCRS> rmap = new HashMap<>();

        for (SingleCRS s : queryCrsComps) {
            for (SingleCRS t : resourceCrsComps) {
                try {
                    CRS.findOperation(s, t, null);
                    if (map.put(s, t) != null) {
                        throw new IllegalArgumentException("Multiple axis can be mapped for " + s);
                    }
                    rmap.put(t, s);
                } catch (FactoryException ex) {
                    //do nothing
                }
            }
        }

        // reduce
        final GridExtent extent = resource.getExtent();
        final long[] low = extent.getLow().getCoordinateValues();
        final long[] high = extent.getHigh().getCoordinateValues();
        for (int resIdx = 0; resIdx < resourceCrsCompsOffset.length; resIdx++) {
            final SingleCRS resourceCrsComp = resourceCrsComps.get(resIdx);
            final SingleCRS queryCrsComp = rmap.get(resourceCrsComp);
            if (queryCrsComp != null) {
                final int qryIdx = queryCrsComps.indexOf(queryCrsComp);

                //intersect with original grid
                final GridGeometry querySlice = reduce(query, queryCrsCompsOffset[qryIdx], queryCrsComp.getCoordinateSystem().getDimension());
                final GridGeometry resourceSlice = reduce(resource, resourceCrsCompsOffset[resIdx], resourceCrsComp.getCoordinateSystem().getDimension());
                final GridExtent intersection = resourceSlice.derive().rounding(mode).subgrid(querySlice).getIntersection();

                for (int i = 0, n = intersection.getDimension(), k = resourceCrsCompsOffset[resIdx]; i < n; i++) {
                    low[k + i] = intersection.getLow(i);
                    high[k + i] = intersection.getHigh(i);
                }
            }
        }

        return resource.derive().resize(new GridExtent(null, low, high, true)).build();
    }
{code}

{code:java}
private static GridGeometry reduce(GridGeometry grid, int offset, int size) {
        final int[] range = new int[size];
        for (int i = 0; i < range.length; i++) {
            range[i] = offset + i;
        }
        return grid.reduce(range);
    }
{code}

{code:java}
private static GridGeometry compound(GridGeometry first, GridGeometry second) throws FactoryException {

        //compound crs
        final CoordinateReferenceSystem compoundCrs = CRS.compound(
                first.getCoordinateReferenceSystem(),
                second.getCoordinateReferenceSystem());

        //compound exents
        final GridExtent secondExtent = second.getExtent();
        GridExtent compoundExtent = first.getExtent();
        for (int i = 0; i < secondExtent.getDimension(); i++) {
            compoundExtent = compoundExtent.insert(
                    compoundExtent.getDimension(),
                    secondExtent.getAxisType(i).orElse(null),
                    secondExtent.getLow(i),
                    secondExtent.getHigh(i),
                    true);
        }

        //compound grid to crs
        final MathTransform compoundTrs = MathTransforms.compound(
                first.getGridToCRS(PixelInCell.CELL_CENTER),
                second.getGridToCRS(PixelInCell.CELL_CENTER));

        //compound grid geometry
        return new GridGeometry(compoundExtent, PixelInCell.CELL_CENTER, compoundTrs, compoundCrs);
    }
{code}

> GridGeometry.subgrid(…) should accept an area of interest with less dimensions than the base grid
> -------------------------------------------------------------------------------------------------
>
>                 Key: SIS-514
>                 URL: https://issues.apache.org/jira/browse/SIS-514
>             Project: Spatial Information Systems
>          Issue Type: Task
>          Components: Coverage
>    Affects Versions: 1.0
>            Reporter: Martin Desruisseaux
>            Assignee: Martin Desruisseaux
>            Priority: Major
>
> If {{grid}} is a 4-dimensional grid geometry, for example with the following content:
> {noformat}
> GridGeometry
>   ├─Grid extent
>   │   ├─Column:   [0 … 4499] (4500 cellules)
>   │   ├─Row:      [0 … 4250] (4251 cellules)
>   │   ├─Vertical: [0 …   39]   (40 cellules)
>   │   └─Time:     [0 …    0]    (1 cellules)
>   ├─Geographic extent
>   │   ├─Lower:  80°00′00″S  180°00′00″W  2021-03-02T00:00:00Z
>   │   └─Upper:  90°00′00″N  180°00′00″E
> {noformat}
> and if {{areaOfInterest}} is a two-dimensional grid, for example with the following content:
> {noformat}
> GridGeometry
>   ├─Grid extent
>   │   ├─Dimension 0: [0 … 255] (256 cellules)
>   │   └─Dimension 1: [0 … 255] (256 cellules)
>   ├─Geographic extent
>   │   ├─Lower :  01°23′41″N  07°52′24″W
>   │   └─Upper :  46°29′05″N  37°12′59″E
> {noformat}
> A call to {{grid.derive().subgrid(areaOfInterest)}} causes an {{OperationNotFoundException}} to be thrown from the {{org.apache.sis.coverage.grid.CoordinateOperationFinder.changeOfCRS}} method because a transform from 2D to 4D can not be found.
> The {{subgrid}} operation should be more tolerant, by leaving untouched the grid in the dimensions not covered by {{areaOfInterest}}.



--
This message was sent by Atlassian Jira
(v8.3.4#803005)