WCPS_Output of CrsTransform() should update all relate metadata from input crs

[Bang Pham Huu]

From Vlad's finding the output of CrsTransform() still in original CRS (i.e: related metadata (geo-axes, crs, bounding boxes, resolutions,..)).

The output of CrsTransform() should update coverage's related metadata to target CRS then user can do subsets on this output from input CRS, e.g:

for $t in (ecmwf_2T) return encode(
  crsTransform(
    $t[ansi("2018-01-31T18:00:00.000Z")], 
    { Lat:"http://www.opengis.net/def/crs/EPSG/0/32632", Long:"http://www.opengis.net/def/crs/EPSG/0/32632"}, {}
  )[E(275985:582315), N(5275985:5582315)], 
"application/gml+xml")

with error

<ows:ExceptionText>Invalid axis name 'E'.</ows:ExceptionText>

[Bang Pham Huu]

This feature will have a problem to get the correct grid bound for geo X-Y axes and update to coverage WCPS metadata, I don't have another idea than have to send rasql query to get the grid bound from Rasql, e.g:

SELECT sdom(
 project( c[1,0:100,0:231], "-40.5,25,75.5,75.5", "EPSG:4326", "EPSG:3785" ), ) FROM test_eobstest_2018_05_09_08_33_46_1355 AS c

Why is the grid bounds not as same as the original ones? because depending on the target CRS, gdal will reproject the 2D domain to different outputs, example:

• source CRS is EPSG:4326, target CRS is EPSG:3857

  $ gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 mean_summer_airtemp.tif output.tif
  Creating output file that is `835P x 770L`.
  

• source CRS is EPSG:4326, target CRS is EPSG:32652

  gdalwarp -s_srs EPSG:4326 -t_srs EPSG:32652 mean_summer_airtemp.tif output1.tif
  Creating output file that is `966P x 831`
  

[Dimitar Misev]

How does subsettingCrs work? It has to determine the grid coordinates for a target CRS. Isn't it what you want to do here?

[Bang Pham Huu]

Replying to dmisev:

How does subsettingCrs work? It has to determine the grid coordinates for a target CRS. Isn't it what you want to do here?

subsettingCRS means translation the input subsets from target CRS to source CRS (i.e: original grid bounds of coverage don't change like with crsTransform()).

[Dimitar Misev]

Ok, maybe check in the GDAL documentation if the reprojection can be done on metadata level (in petascope via GDAL JNI).

[Bang Pham Huu]

Replying to dmisev:

Ok, maybe check in the GDAL documentation if the reprojection can be done on metadata level (in petascope via GDAL JNI).

There is ​`gdal.ReprojectImage()` which can be used in Java GDAL JNI. It still needs input dataset though (so a dummy 2D dataset (e.g: with only 0 value) needed to be created before this function) then can get the domains of output.

[Dimitar Misev]

That won't work, the input has to match exactly the width/height.

[Bang Pham Huu]

Replying to dmisev:

That won't work, the input has to match exactly the width/height.

What do you mean? we know the grid bounds of the original one and we create the dummy with this size as input for gdal ReprojectionImage().

[Vlad Merticariu]

We don't want to create a dummy dataset, this can be potentially large. You can use other libraries to determine the new bounds (maybe geotools).

[Dimitar Misev]

Right, but this is still very inefficient. You'll need to create a 1GB array if the data is 10,000 x 10,000 and gdal will work to reproject it, just to get the resulting width and height.

[Bang Pham Huu]

Replying to dmisev:

Right, but this is still very inefficient. You'll need to create a 1GB array if the data is 10,000 x 10,000 and gdal will work to reproject it, just to get the resulting width and height.

I know but no choice as GDAL wants to have input. This is resampling data problem so it must needs neighbor pixels to interpolate for example, so just naive input bbox cannot work.

[Dimitar Misev]

I doubt that's the case, because you likely need to initialize the output array with the correct result width/height already before doing the resampling.

[Vlad Merticariu]

Check out proj4 and geotools for a start:

​http://neondataskills.org/GIS-spatial-data/Working-With-Rasters/

​http://docs.geotools.org/latest/tutorials/geometry/geometrycrs.html

I see that geotools can reprojects vectors, so if you define a polygon that encompasses the coverage and reproject that one you can find out the new extents. Probably not a good solution, but it shows that it can be done without touching the data values.

[Dimitar Misev]

From: ​http://www.gdal.org/gdalwarper_8h.html#ab5a8723d68786e7554f1ad4c0a6fa8d3

The ​GDALSuggestedWarpOutput() function is used to determine the bounds and resolution of the output virtual file which should be large enough to include all the input image

[Dimitar Misev]

But still needs the input dataset it seems

[Dimitar Misev]

I looked at the code of this method and it doesn't seem to touch the data actually; it just uses the width/height from hSrcDS.

[Bang Pham Huu]

Replying to vmerticariu:

Check out proj4 and geotools for a start:

​http://neondataskills.org/GIS-spatial-data/Working-With-Rasters/

​http://docs.geotools.org/latest/tutorials/geometry/geometrycrs.html

I see that geotools can reprojects vectors, so if you define a polygon that encompasses the coverage and reproject that one you can find out the new extents. Probably not a good solution, but it shows that it can be done without touching the data values.

It doesn't work, the way it converts for each vertex in polygon is as same as we've done in clip extension. You don't know the pixel sizes after all.

The real problem here is about resampling data and no way to do it as it requires the input dataset (GDAL, geotools or whatever).

[Vlad Merticariu]

Are you saying that the pixel value matters in reprojection?

[Bang Pham Huu]

Replying to vmerticariu:

Are you saying that the pixel value matters in reprojection?

https://image.slidesharecdn.com/gisconcepts3-091126070346-phpapp01/95/gis-concepts-35-53-728.jpg?cb=1259219087

check it out, different resampling method will use different equations.

[Vlad Merticariu]

Yes, resampling is necessary for finding out the new pixel values, not for finding out the extents or the pixel sizes.

[Bang Pham Huu]

Replying to vmerticariu:

Yes, resampling is necessary for finding out the new pixel values, not for finding out the extents or the pixel sizes.

it must! resample before it can find the extent as it cannot imagine without calculation. Check also with projection

​http://www.mdpi.com/ijgi/ijgi-06-00092/article_deploy/html/images/ijgi-06-00092-g001.png

[Vlad Merticariu]

In none of these examples I see the value of the pixels being used. The shape of the input raster you know in petascope.

I'm maybe missing something, but where exactly is the value of the pixels used to determine the shape of the output?

If the value is not used anywhere, then you don't need a raster, you only need it's shape which is metadata.

[Bang Pham Huu]

Replying to vmerticariu:

In none of these examples I see the value of the pixels being used. The shape of the input raster you know in petascope.

I'm maybe missing something, but where exactly is the value of the pixels used to determine the shape of the output?

If the value is not used anywhere, then you don't need a raster, you only need it's shape which is metadata.

check it again ​http://www.mdpi.com/ijgi/ijgi-06-00092/article_deploy/html/images/ijgi-06-00092-g001.png

see how the output pixels is created based on the CRS projection.

[Vlad Merticariu]

Again, where is the value of the pixel (e.g. 255c) used in that? In the example only the location of the pixel is used.

[Bang Pham Huu]

Replying to vmerticariu:

Again, where is the value of the pixel (e.g. 255c) used in that? In the example only the location of the pixel is used.

so in the output example, what is the value of target pixels you think it will fill? Also, without filling pixel by pixel for the output image, how the output extent should look like?

[Vlad Merticariu]

The value of the target pixel is the value of the input pixel, but it does not matter what this value is.

The problem itself is the following: given the grid domain, geo domain, resolution and crs of an input raster and a target crs into which the raster is to be projected, find out the shape (grid domain, geo domain, resolution) of the output raster.

Your task is to find solutions to this problem.

I agree that one solution is to simply take every pixel and reproject it, then take the min / max, but this is not efficient. I don't believe that this is the only solution, at least in some cases it could be avoided. As Dimitar suggested above, gdal does this incrementally. Check the gdal implementation below:

​https://github.com/OSGeo/gdal/blob/106c8288e7a05f4efc1a588c5a3b2da7ec52d915/gdal/alg/gdaltransformer.cpp#L354

If you have problems understanding the syntax, let me know (line(s) that you don't understand).
Also check with Brennan, he might be able to give you some pointers.

[bbell]

Replying to vmerticariu…

So, I see a few variables, grid sdom s, geo sdom g, some mystery function T(g) = g', which transforms the geo sdom to some g', and the problem is that we want to know s', but our only add'l information is the transformation from g to s, let's say R(g) = s.

If we had R': g' → s', we could figure this out straight away since R' T R^-1(s) = s'

So, I read through the documentation of gdaltransform a bit, because I wanted to know how R' is determined.

 * Then a resolution is computed with the intent that the length of the
 * distance from the top left corner of the output imagery to the bottom right
 * corner would represent the same number of pixels as in the source image.
 * Note that if the image is somewhat rotated the diagonal taken isn't of the
 * whole output bounding rectangle, but instead of the locations where the
 * top/left and bottom/right corners transform.  The output pixel size is
 * always square.  This is intended to approximately preserve the resolution
 * of the input data in the output file.

Since R' is just a resolution scaling, all that you need to compute is T(x₀) and T(x_f) where g = [x₀:x_f].

[Bang Pham Huu]

Thanks for all your's helps, the ticket was solved with using gdal-java version 1.10.1