Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Get data value from a point #112

Open
jratike80 opened this issue Aug 11, 2020 · 17 comments
Open

Get data value from a point #112

jratike80 opened this issue Aug 11, 2020 · 17 comments
Assignees

Comments

@jratike80
Copy link

Users might want to query a coverage for getting a data value from a single point. Typical use case is to check the height at given location from DEM but probably queries with more dimensions "temperature at location x,y,z at time …" make sense as well.

WCS 2.0 defines two possibilities for getting a coverage from a single point.

  1. Use subsets with lower limit = upper limit (trim to a point)
    &SUBSET=Lat(20,20)&SUBSET=Long(30,30)

  2. Use subsets with slice points (slice to a point)
    &SUBSET=Lat(20)&SUBSET=Long(30)

The difference between the trimming and slicing is that trimming does not alter the dimensionality but slicing removes dimensions. If the coverage is 2D raster DEM in case 1) the result can be 1 by 1 pixel sized GeoTIFF. In case 2) the result must use format like GML because rasters have always 2 dimensions but the result must have zero.

Behavior in existing WCS 2.0 implementations vary. GeoServer sends 1x1 sized GeoTIFF with slice type of query but trim to point yields an error "Empty intersection after subsetting"
https://osgeo-org.atlassian.net/browse/GEOS-9553. MapServer sends 1x1 sized GeoTIFF with trim to point query but it does not accept slicing query. Because MapServer is limited to send coverages only as rasters it cannot use for example GML for 0-dimensional coverages.

To my knowledge it is much more common to use WMS GetFeatureInfo for pixel value queries than WCS. However, GetFeatureInfo does not suit well for other than 8 bit rasters, and if raster is styled for visualisation, some servers send the pixel value of the native data but others the pixel value after styling is applied. Thus GFI is a poor workaround and users deserve better.

I would like to see a recommendation about how to support point queries in OACov. From the client application side I believe that slicing together with some very simple format for rangeset would be convenient.

Unfortunately point queries are not as trivial as they may feel. In the next image query by point 1 that hits a data point should not give any troubles. Point 2 is also easy and those who think that coverages are pixels can select the intersected pixel, and folks who prefer points can take the nearest point. The result is the same.

Queries by points 3 and 4 are not as easy. Should the result be an average of 2 or 4 nearest points, respectively (requires interpolation), or simply the value of some nearest neigbor? Which nearest? How to guarantee that different servers give uniform results? There are also coverages which are certainly not continuos and only query by point 1 is reasonable.

point_query

@jerstlouis
Copy link
Member

jerstlouis commented Aug 11, 2020

0D support (slice to point) could return single value range set, which in application/json media-type could be either a number for a single band e.g.

403.51

or an object for multiple bands:

{ "red" : 0.434, "green" : 0.123, "blue" : 0.53 }

Trim to point could return 1x1 image.

@joanma747
Copy link

Or could be a CIS document that specifies that the domainset is a point in space (and the coordinates, and the rangetype clarifies the meaning of "red" "green" and "blue" and rangeset is something similar to the object above.

@cmheazel
Copy link
Contributor

This content has been captured in the API-Coverages Users Guide. Feel free to expand and refine.

@ghobona ghobona transferred this issue from opengeospatial/OGC-API-Sprint-August-2020 Nov 24, 2020
@chris-little
Copy link

Should the result be an average of 2 or 4 nearest points, respectively (requires interpolation), or simply the value of some nearest neigbor? Which nearest? How to guarantee that different servers give uniform results? There are also coverages which are certainly not continuos and only query by point 1 is reasonable.

@jratike80 The optimal interpolation is certainly determined by the result/variable/parameter of interest. E.g. surface pressure versus rainfall. Interpolation is reasonable for pressure on small scales, whereas choosing the nearest point value is best for most rainfall coverages.

@percivall
Copy link
Collaborator

percivall commented Dec 2, 2020 via email

@chris-little
Copy link

Abstract spec topic 6
Annex C Interpolation methods

Spotted a typo in third paragraph of C.6 - 'p' should be 'q' !

@pebau
Copy link
Contributor

pebau commented Dec 6, 2020

@chris-little : thanks for the hint, Chris - just trying to verify (important for our ISO 19123-1 work), but in C.6 Bilinear interpolation I cannot find any p nor q. Can you help me spotting it?

@chris-little
Copy link

chris-little commented Dec 9, 2020

Typo by me. It was in C.9. p should be q. Probably doesn't affect anyone's code.

C.9Lost area interpolation
Lost area interpolation extends a discrete point coverage to a continuous function, f, defined on the convex hull of the domain of the point coverage.Let {p1, p2, ..., pn} be the domain of the point coverage, and let {T1, T2, ..., Tn} be the Thiessen polygons generated by the set {p1, p2, ..., pn}.Suppose it is desired to calculate f(q), where qis a direct position in the convex hull of {p1, p2, ..., pn}. Begin by forming the Thiessen polygons generated by {p1,p2, ..., pn}; then add pto the set {p1, p2, ..., pn}, and form the Thiessen polygons for the set of n+1 points: {p1, p2, ..., pn, q}. The two sets of polygons are identical, except that each of the polygons coterminous with the polygon containing q―loses area‖ to the new polygon containing q.

@jerstlouis jerstlouis added this to Backlog in Core release via automation Feb 3, 2021
@jerstlouis
Copy link
Member

jerstlouis commented Feb 2, 2022

SWG 2022-02-02: Suggesting that if slicing to 0 dimension (e.g. specifying a single value for both latitude and longitude for a 2D geospatial coverage):

  • If the negotiated format only supports 2D, then the 1x1 value should still be returned as a 2D coverage,
  • If the negotiated format supports returning a single value, then the response can be 0D
    • NOTE: CIS JSON currently doesn't support returning a 0D coverage, does it?

Note same would apply if user slices to e.g. a single line and negotiated format only supports 2D. CIS JSON in this case would only have one axis in the response.

@jerstlouis jerstlouis moved this from Backlog to In Progress in Core release Feb 2, 2022
@pebau
Copy link
Contributor

pebau commented Feb 2, 2022

@jerstlouis this feels like a solution for one specific situation only, but not a general one:

  • can a 2D only format represent a 1D result?
  • TIFF supports single values, is that 2D or 0D? How should TIFF differentiate between trimming and slicing when both lead to a single value remaining?

Further, isn't that a physical issue = determined by the concrete encoding format?

@pomakis
Copy link

pomakis commented Feb 2, 2022

For what it's worth, CubeWerx's implementation of the "OGC API - Coverage" endpoints does as @jerstlouis suggests when serving coverages as GeoTIFFs. That is, if a 1D slice is requested, we return a 1xn or nx1 GeoTIFF, and if a 0D slice is requested, we return a 1x1 GeoTIFF. That's the only behavior that makes sense to me. And @jerstlouis's generalization of this to any 2D-only format seems reasonable and practical.

The question of whether 1x1 GeoTIFFs are 2D or 0D seems needlessly pedantic, and it goes down a rabbit hole. I mean, wouldn't that technically depend on whether the coverage is PixelIsPoint or PixelIsArea? I think it's reasonable for the specification to state, without getting tangled up in theoretics, that for 2D-only formats such as GeoTIFF, a 1D slice is represented as a single row or column, and a 0D slice is represented as a single cell.

@jratike80
Copy link
Author

jratike80 commented Feb 2, 2022

My point when I created this issue was that both we ourselves at NLS and our customers need to do such queries from our WCS service. Usually with our data the interest is in height at a point or height profile along a path. Most of our users seem to prefer asciigrid over GeoTIFF as a format for points. Measures along a line they probably collect on client side by doing first a bunch of smallish asciigrid or GeoTIFF GetCoverage requests. I suppose a direct 1D response would be appreciated for profiles but our server can output only raster formats.

Would it be better to have an own request type like GetFeatureInfo is in WMS, perhaps not in core but as an extension? That would allow to define separate response formats for queries. WMS layers have also metadata queryably yes/no that might be reasonable as well for coverages, at least when the response to a query by point or by linestring would require either interpolation or selection of the nearest points. @chris-little pointed out that interpolation may suit for some data but nearest value for some other, and I guess that some data may not be suitable for either.

@jerstlouis
Copy link
Member

@pebau

this feels like a solution for one specific situation only, but not a general one

We are trying to make this general enough to handle most encodings, and at least 0D and 1D responses.

can a 2D only format represent a 1D result?

Yes it can by the extra dimension axis having only one possible value along its axis. It is not expressly 1D, but embedded metadata and/or the context can clearly inform how to interpret it as 1D. This is what we are trying to clarify here, perhaps with a clearly specified permission in the specification.

TIFF supports single values, is that 2D or 0D?

TIFF itself as a format is limited to 2D, but this permission would allow to also understand and properly interpret 2D responses as 0D or 1D for a slicing request which would result in a 0D or 1D response if not for the format limitation.

How should TIFF differentiate between trimming and slicing when both lead to a single value remaining?

At the moment it cannot, but a future GeoTIFF version embedding e.g. CIS JSON could make it possible.

Further, isn't that a physical issue = determined by the concrete encoding format?

Yes, but with a permission we can clearly specify how this can still work for any formats limited to 2D.

@jerstlouis
Copy link
Member

@jratike80 I like the pun ;)

The height profile along a path use case is exactly what the EDR trajectory query allows to do. Perhaps a geometry query could also be integrated with the Coverages API at some point. In this case it would make sense to return this as 1D, in whichever negotiated format (be it GeoTIFF, or ASCIIGrid, or CIS JSON or anything else).

As for interpolating or nearest, as I mentioned today I feel this should always follow the nature of the coverage range type, where each field in the data record either makes sense to interpolate or not, and servers should behave in accordance with that.

@cnreediii
Copy link
Contributor

For the CDB analytics use cases defined for version 2.0 (as per the 2020 CDB Sprint), getting a value at a point in a coverage is required for numerous analytic operations such as cross sections and profiles.

@jerstlouis
Copy link
Member

SWG 2022-02-16: We'll add the permission to return a response in more dimension than the "sliced" results for formats with a minimum number of dimensions (e.g. a 1x1 GeoTIFF for requesting a single point).

For the trajectory use case, we will consider an extension for specifying geometry (#52) and possibly take some inspiration from the EDR API trajectory queries. How should the response look like in CoverageJSON / CIS JSON / netCDF? What is the domain of the response? Will the resulting coverage be a "point cloud" with coordinates for each point, or will it simply use an index CRS corresponding to the points in the geometry of the request?

@jerstlouis
Copy link
Member

We should consider defining a GeoJSON conformance class (#166) which would easily and unambiguously support this use case of querying a single value (Point) or for relatively small point clouds (MultiPoint).
The properties would correspond to the fields (range), and the point geometry would provide the domain.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
Core release
  
Agreed; to be applied
Development

No branches or pull requests

10 participants