Recognize if an image is geo-referenced or not with GDAL

[Bang Pham Huu]

When i'm testing Rasdaman by using script test.sh in test_wcps, it will import rgb file by rasql query in file petascope.sh, but I also using WCST import this file with map_mosaic recipe.

It seems that WCST import tool work good when importing non-georeferened image like RBG because it will take the origin is the top left (0, 0).

When using rasql query it will choose the origin is the lower left (0, 0).

So I guess something not true with rasql importing. See the attachment for more detail.

[Bang Pham Huu]

Difference between using wcst-import (rgb1) and rasql (rgb)

[Dimitar Misev]

So which one is correct? wcst_import.sh or rasql?

[Bang Pham Huu]

rasql is not correct, so that is why I add cc to you, Dimitar.

[Bang Pham Huu]

I can, Dimitar, from my point of view, the origin should be on the top left and user should query from this origin (top to bototm) instead of rasql (origin in lower left) (bottom to top).

We could discuss more on Wednesday for more detail, thanks.

[Alex Dumitru]

Both methods are correct (for non-geo-referenced images), the choice of choosing the origin in the OGC services is left to the user. Btw, this has nothing to do with rasql, it only has to do with the origin metadata in petascopedb.

WCSTImport by default sets the origin according to the MetadataProvider used(which in all current recipes is GDAL which considers the origin always as top left corner). I doubt that there is any reference in the WCPS standard of what the origin should be, but someone should check it just in case it's defined there.

[Dimitar Misev]

So this only affects non-geo-referenced data right, and it's because psql statements in the systemtest set the origin differently from wcst_import?

I'd say it's safe to update the oracles then, as long as they are really wrong because of the different origin only.

[Bang Pham Huu]

I will add a details with origins problem below. From my view, I think there should be separation between non-georeferenced (using top-left origin is lowest point) and geo-referenced images (using low-left origin is lowest point).

But as Alex mentioned, there could be configured for choosing the origin is a good idea if user don't familiar when they using image with origin is top-left corner.

In test oracle, mr also has this problem. Please make sure that test oracle use the "origin configuration" or not, I don't know how to check it.

This example here is a PNG rgb width = 400 (x), height = 344 (y)

1 is using rasql + metadata in petascopedb:

 <origin>
          <Point gml:id="rgb-origin" srsName="http://localhost:8080/def/crs/OGC/0/Index2D" axisLabels="i j" uomLabels="GridSpacing GridSpacing" srsDimension="2">
            <pos>0 344</pos>
          </Point>
        </origin>

2 is using wcst-import

It seems the origin is (1, 1) should be (0, 0)

<origin>
          <Point gml:id="rgb1-origin" srsName="http://localhost:8080/def/crs/OGC/0/Index2D" axisLabels="i j" uomLabels="GridSpacing GridSpacing" srsDimension="2">
            <pos>1 1</pos>
          </Point>

[Dimitar Misev]

geo-referenced images (using low-left origin is lowest point)

why?

But as Alex mentioned, there could be configured for choosing the origin is a good idea if user don't familiar when they using image with origin is top-left corner.

where did Alex mention this?

It seems the origin is (1, 1) should be (0, 0)

yes indeed, I'd expect (0, 0).

[Alex Dumitru]

@Bang, I notice that in your patch you check if a file is not georeferenced with this condition geotransform[0] == geotransform[2] == 0. geotransform[2] will almost always be 0 in gdal, as it only applies for files with strange CRSs and warped axes and geotransform[0] is 0 in many cases (e.g. when the image starts at 0 degrees Longitude). Am I missing something?

@Dimitar Do we have georeferenced images in systemtests as well? It would be nice to test all the recipes that we offer in the systemtests as well. I can provide some data for this if we don't already have some.

[Dimitar Misev]

mean_summer_airtemp should be geo-referenced, it would be nice to have some more such data though.

[Bang Pham Huu]

Replying to dmisev:

geo-referenced images (using low-left origin is lowest point)

why?

it you latitude and longtitude value so (0, 0) is the center and that is why low-left is the lowest point (it is opposite with non-georeference image when (0, 0) is in the top-left).

But as Alex mentioned, there could be configured for choosing the origin is a good idea if user don't familiar when they using image with origin is top-left corner.

where did Alex mention this?

He said this "the choice of choosing the origin in the OGC services is left to the user", and I think he want to mention that user can config when using wcst-import to set the origin of image (for example non-georeference image will have origin in bottom-left instead of top-left) ———> so this is needed to separate 2 different cases (I still think this should not be chosen by user but from original definitions of origin with non-georeferenced and geo-referenced image).

It seems the origin is (1, 1) should be (0, 0)

yes indeed, I'd expect (0, 0).

Yes it should be (0, 0), because

i (measured in GridSpacing ) with domain extent from 1 to 399
j (measured in GridSpacing ) with domain extent from 1 to 343

looks like user will lost a column of pixels on the left most side.

[Bang Pham Huu]

Replying to mdumitru:

@Bang, I notice that in your patch you check if a file is not georeferenced with this condition geotransform[0] == geotransform[2] == 0. geotransform[2] will almost always be 0 in gdal, as it only applies for files with strange CRSs and warped axes and geotransform[0] is 0 in many cases (e.g. when the image starts at 0 degrees Longitude). Am I missing something?

I need to add this condition because from ​http://rasdaman.org/ticket/1018, you will see that when using wcst-import non-georeferenced files like RGB, mr,…the j axis will be wrong order like (344, 0) instead of (0, 344). I don't mention about the georeferenced images which your code was working correctly. Please note that also, your code will have problem with the float origin of non-georeference image (​http://rasdaman.org/ticket/941#comment:20)
For example: mean_summer_airtemp is this kind of image it is non-georeferenced but the origin is not basically (0,0).

Upper Left  ( 111.9750000,  -8.9750000) 
Lower Left  ( 111.9750000, -44.5250000) 
Upper Right ( 156.2750000,  -8.9750000) 
Lower Right ( 156.2750000, -44.5250000) 
Center      ( 134.1250000, -26.7500000

@Dimitar Do we have georeferenced images in systemtests as well? It would be nice to test all the recipes that we offer in the systemtests as well. I can provide some data for this if we don't already have some.

I can answer this question, yes we have files in irr_cube_2 that has georeferenced images. We also have 2 kind of strange files (that is non-georeferenced but origin is not (0, 0)) is eobs.nc (NetCDF File) and (mean_summer_airtemp.tif). And if we don't set reference to these images, using wcst-import tool will raise error ​http://rasdaman.org/ticket/941#comment:19

So, I will leave a question for both of you that there a case non-georeferenced images will have the origin is not (0, 0) and we should know how to solve this problem (it is easy to set geo-referenced for this kind of images but if user want to store the original images but not relate with geography and don't care with CRS, it will noot good for them).

[Alex Dumitru]

@Bang: I still don't get why the geotransform[0] == geotransform[2] == 0 actually means that the image is not georeferenced? How did you figure out this is the case? To me it seems this condition is valid for some geo-referenced images as well.

[Bang Pham Huu]

Hi Alex,

It is very basic knowledge of Computer Graphic when image (generally without geo-referenced) is used the top-left as origin (lowest point). See an example here from GDAL

Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  720.0)
Upper Right (  720.0,    0.0)
Lower Right (  720.0,  720.0)
Center      (  360.0,  360.0)

and the answer for why geotransform[0] == geotransform[2] == 0 is non-georeferenced image is

# My comment for better understanding
 GDAL Geotransform:
          + 0: top left x
          + 1: west-east pixel resolution (x axis)
          + 2: 0 value
          + 3: top left y
          + 4: 0 value
          + 5: north-south pixel resolution (y axis)

And for geo-referenced image when you use Cartesian coordinates with X axis Longtitude and Y axis Latitude, the lowest point (origin) will change to bottom-left, see below

Upper Left  (-110.0000000,  80.0000000) (110d 0' 0.00"W, 80d 0' 0.00"N)
Lower Left  (-110.0000000,  20.0000000) (110d 0' 0.00"W, 20d 0' 0.00"N)
Upper Right ( -60.0000000,  80.0000000) ( 60d 0' 0.00"W, 80d 0' 0.00"N)
Lower Right ( -60.0000000,  20.0000000) ( 60d 0' 0.00"W, 20d 0' 0.00"N)

This is why your code not work and need my code for the case of "non-georeference image". It is best if you try importing "non-georeferenced image" by wcst-import and see exception.

[Alex Dumitru]

@Bang: You claim geotransform[2] is always 0, correct? You also claim that if geotransform[0] == 0 then the image is not geo-referenced. Wouldn't a image georeferenced with these coordinates be a non-geo-referenced image applying your rule?

Upper Left  (0,  80.0000000) (...)
Lower Left  (0,  20.0000000) (...)
Upper Right (20.0000000,  80.0000000) ( ... )
Lower Right (20.0000000,  20.0000000) ( ... )

The origin is not the lowest or highest point, it's just some point that somebody chooses. The offset vector then tells you how to navigate to other points. As far as I know, the origin of the Petascope Coverages that WCSTImport ingests is acutally in the upper left corner and not in the bottom left.

I know that wcst_import does not work correctly with non-geo-referenced images, I'm just trying to understand how to fix it.

[Bang Pham Huu]

yes, you're quite right with my example (actually if georeferenced image with this kind of coordinates will be non-georeferenced image).

But if we could get the coordinate of top-left and bottom-left, we could check correctly, it is not difficult I think. From my previous comment, I suggest that (for example X is width and Y is height). Top-left is X1,Y1 and Lower-left is X2, Y2.

so if (Y1<Y2) then we are importing non-georeferenced image.
and if (Y1 > Y2) then we are importing geo-referenced image.

Sorry for making use confuse last time with wrong condition but my idea still the same from the origin of non-georeferenced and georeferenced. Now I don't think we could prove it is wrong.

[Peter Baumann]

I see your point, Bang, but I wonder whether we can deduce georeferencing - is "0 at top" indicative, or is it just a convention used mostly, but not always? Also consider that we go for Planetary cases → maybe asking Angelo would be good.
So the goal for me seems: establish a common rul on how to determine whether a coverage is georeferenced. In case of a complete coverage we have CRS info which tells us, so this discussion seems to focus specifically on the case that CRS info is not available.
I understand that >1>Y2 indicates some georeferencing, but does this help us if we don't know a CRS?
Should it turn out that this is not decidable (which would make sense to me) we just would need to resort to a GridCoverage with IndexCRS only.
Curious about the outcome.

[Alex Dumitru]

I would argue, that at least from the point of WCSTImport, which we are discussing now, the only indicative that a coverage is georeferenced is having a CRS and that CRS is different than /def/crs/OGC/0/Index{SomeValue}D
Every other method that is based on gdal geotransform has plenty of edge cases when a georeferenced image would be identified as non-georeferenced.

With that in mind, I still consider that we are not seeing the full picture here. The gdal documentation says (​http://www.gdal.org/gdal_datamodel.html#gdal_datamodel_dataset_gtm) that this formula should give you the georeferenced x,y coords based on the raster x,y coords:

 maxx = minx + geo_transform[1] * self.gdal_dataset.RasterXSize
 miny = maxy + geo_transform[5] * self.gdal_dataset.RasterYSize

Why does it fail for non-georeferenced images? We should ask on ​http://gis.stackexchange.com/ or on the gdal mailing list.

[Peter Baumann]

I'm lost: if it is non-georeferenced, why do you want to apply the formula? Isn't in this case the pixel extent to be used directly?

[Alex Dumitru]

Well, that's the problem, that GDAL offers no method of determining if the dataset is georeferenced or not and based on the documentation the formula should work for all cases.
We could check the CRS and consider it is not georeferenced if there is none. But that doesn't cover all our cases. If the image has a CRS and it is Index2D it still behaves as a non-georeferenced image so we have to check for this as well. And what happens when someone else defines a similar crs? Do we check for that as well?

I would rather have a general method that works in all cases, and I find it quite strange that GDAL returns the top left corner in (geotransform[0], geotransform[3]) for geo-referenced images and the bottom left corner in the same variables for non-georeferenced images and this is not documented anywhere. That's why I am thinking that we might be missing something here and don't want to rush and implement something that seems to work at a first glance but we don't have any idea why it does.

Either way, Bang's proposed formula for checking if the image is georeferenced (geotransform[0] == geotransform[2] == 0) is not correct and I showed in comment:15 a counter example to it.

[Alex Dumitru]

Replying to bphamhuu:

yes, you're quite right with my example (actually if georeferenced image with this kind of coordinates will be non-georeferenced image).

But if we could get the coordinate of top-left and bottom-left, we could check correctly, it is not difficult I think. From my previous comment, I suggest that (for example X is width and Y is height). Top-left is X1,Y1 and Lower-left is X2, Y2.

so if (Y1<Y2) then we are importing non-georeferenced image.
and if (Y1 > Y2) then we are importing geo-referenced image.

This doesn't work either, now you are assuming that we always have a negative offset vector (geotransform[5]) for geo-referenced images which I don't think is always the case (I can define a coverage that does not fulfill this requirement easily). It might be the case that GDAL always transform the Y-axis offset vector into a negative one and then your algorithm would be correct, but I see no hint of such a thing being true in the documentation (if you point me to it, I'd be happy, so we can finally get a solution).

Sorry for making use confuse last time with wrong condition but my idea still the same from the origin of non-georeferenced and georeferenced. Now I don't think we could prove it is wrong.

[Bang Pham Huu]

Replying to mdumitru:

I would rather have a general method that works in all cases, and I find it quite strange that GDAL returns the top left corner in (geotransform[0], geotransform[3]) for geo-referenced images and the bottom left corner in the same variables for non-georeferenced images and this is not documented anywhere. That's why I am thinking that we might be missing something here and don't want to rush and implement something that seems to work at a first glance but we don't have any idea why it does.

​http://www.gdal.org/gdal_tutorial.html here is where it is documented (Getting Dataset Information)

This doesn't work either, now you are assuming that we always have a negative offset vector (geotransform[5]) for geo-referenced images which I don't think is always the case (I can define a coverage that does not fulfill this requirement easily). It might be the case that GDAL always transform the Y-axis offset vector into a negative one and then your algorithm would be correct, but I see no hint of such a thing being true in the documentation (if you point me to it, I'd be happy, so we can finally get a solution).

The best explanation is from examples, so I will describe problem again and please give me some exception with detail also (sorry if it could take your time for creating image if you think it is not right).

From the problem when I imported non-georeferenced image like rgb, mr (​http://rasdaman.org/ticket/1018) it will have problem with subset j (subset2=j(344.0,0.0)) it should be (0.0, 344.0). And I claim that we should enhanced your code to correct this problem with non-georefrenced image.

I would point out my latest formula has an exception case (please note that I do not assume we are using "negative or positive" offset vector). Just think we have 4 points named A (upper left) and B (lower left) C (upper right) and D (lower right) respectively.

In here I'm using A(x1, y1) - Top-left and B(x2, y2) - Lower-left. And my formula is:

if y1 < y2:
   cout << "non-georeferenced"
else
   cout << "geo-referenced"

+ Non-georeferenced images like rgb, png with top-left (0,0) then y1 = 0 will < y2 = 400 → true.
+ Geo-referenced images like geotiff with top-left (here with negative or positive Y) like:

Y1=0, Y2 = -500 or Y1 = 100, Y2 = 0 or Y1 = -200, Y2 = -300 → true

But we have an exception with non-georeferenced image (and here that is beyond my knowledge - we should notice this case very carefully), because this will prove my formula is wrong. It is file NetCDF ebos.nc in system test of Rasdaman. This is file non-georeference like but the origin is the same as geo-referenced image.

gdalinfo '/home/rasdaman/eobs.nc' 
Driver: netCDF/Network Common Data Format
Files: /home/rasdaman/eobs.nc
Size is 232, 101
Coordinate System is `'
Origin = (-40.500000000000000,75.500000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)

Corner Coordinates:
Upper Left  ( -40.5000000,  75.5000000) 
Lower Left  ( -40.5000000,  25.0000000) 
Upper Right (  75.5000000,  75.5000000) 
Lower Right (  75.5000000,  25.0000000) 
Center      (  17.5000000,  50.2500000) 

In conclusion, I agree that we could not conclude which is geo-referenced or non-georeferenced image with just 4 coordinates. But to fix ticket 1018 to import non-georeferenced image with origin (0,0) in top-left will work with my very first formula but now change from if ( X1 == Y1 == 0) to:

if Y1 < Y2  (0, 0) (top - left non-georeferenced images)
   minyy = Y1
   maxyy = Y2
else  (other cases of non-geoferenced or geo-referenced images)
   minyy = Y2
   maxyy = Y1

I also want to suggest that this should not be overcomplicated by finding out which image is non-georefenced or geo-referenced (because before checking ebos.nc - I always think non-georeferenced image has top-left is (0,0)). It should be concluded that in which case we need to swap minyy and maxyy and that is done.

If we really need to find out what is non-georeferenced or georeferenced, it should be complicated with many types of data that could embedded or non-embedded this information ​http://gis.stackexchange.com/questions/18672/how-to-determine-if-a-tiff-is-georeferenced-or-not

That is all from my points of view and I'm waiting for your decision.

[Bang Pham Huu]

I aslo want to add with this ticket could not solved because when query "rgb imported by wcst-import" like this "for c in rgb( c [ i(0:400), j(0:100) ]" it will return the top-head with origin (0,0) top-left but if query imported by rasql it will return the bottom-tail with origin (0,0) lower-left.

So we should have an identical case in there and I think it is best that image imported by rasql will use origin in top-left (only with this case). Or basically saying, just choose the lower-point of left-side is the origin for logical understanding.

[Peter Baumann]

1 - having an IndexCRS is the formal way of saying "non-referenced". So whenever a coverage is decided to be non-referenced we would pick an IndexCRS to express this.
2 - FWIW, I am sure Bang is aware of it: with all CRSs (including IndexCRS) we can have a negative axis direction to accommodate all possible cases.
3 - Most important: IMHO we do not need to make a forced decision - a coverage must be sufficiently complete. Interpreting incomplete information (such as from a binary file only, without GMLCOV XML) we should not do if there is any slight chance of ambiguity. Just map it to a non-referenced coverage, ie: IndexCRS. Yes, we might lose some information, but it forces users to be comprehensive if they want. If no geo CRS is indicated we anyway cannot make it a referenced coverage at all → No CRS identifier, no georeference

[Dimitar Misev]

So in conclusion, we can simply check if IndexND CRS is used, it seems like that would be sufficient for now.. until we figure out some other non-georeferenced CRSs.

[Bang Pham Huu]

I agree with Dimitar with that conclusion. However, it looks like we has been moved outside of this ticket and ticket ​http://rasdaman.org/ticket/1018 due to my think that this because we must separate non-georeferenced and georeferenced images to solve ticket 1018. As you can see in my comment ​http://rasdaman.org/ticket/1021#comment:22 I've pointed out even non-georeference image could have coordinates as geo-referenced images.

So I hope that everything should not be mixed up. The solution to know which is non-georeferenced or georeferenced image is not related to both tickets. This could be good idea if in some case we need to know to handle but not in this case.

I want to confirm again, that ticket relate to use wcst_import to import file RGB and rasql + metadata petascope in petascope.sh (test_wcps) has different result. I think it should be considered carefully because 1 query could not return different results.

Please see attachment below and reconsider.

​http://rasdaman.org/attachment/ticket/1021/abcd.jpg

[Peter Baumann]

FWIW, there is no other non-referenced CRS besides IndexCRSs, as these have been established sepcifically for expressing this situation.

[Peter Baumann]

Bang, you are right: we need to come to a clear decision on the axis direction, in particular for the IndexCRS. This is not yet clear from standardization, and any choice we make will have pros and cons. My naive idea is that an IndexCRS, reflecting a Cartesian CRS, has origin _bottom_ left.

As I am not sure myself I have asked OGC lists.

[Dimitar Misev]

Replying to pbaumann:

Bang, you are right: we need to come to a clear decision on the axis direction, in particular for the IndexCRS. This is not yet clear from standardization, and any choice we make will have pros and cons. My naive idea is that an IndexCRS, reflecting a Cartesian CRS, has origin _bottom_ left.

As I am not sure myself I have asked OGC lists.

Let me add Piero as the creator of the IndexND CRSs, I guess he would have a quick answer.

[Dimitar Misev]

In the CRS defs, axis direction is ​positive and origin is the ​lower bound of the axis

[Alex Dumitru]

The link you provided does not say that the origin is in lower bound, or maybe I'm missing the point somehow. The definition says that the *origin of the coverage* can be anywhere, but the *origin of the geopixel* is in the lower corner of the geopixel as opposed to others where it is in the middle of the geopixel.

As far as I know, nothing in the CRS ever dictates the origin of the coverage. The origin can be wherever you want it to be. The only responsibility that you have as a coverage creator is that origin + offset_vector * number_of_raster_pixels = extent

[Dimitar Misev]

Ah right, in the description I misinterpreted "lower bound of the cell" to be about the origin of the coverage.

[Dimitar Misev]

Response from Meteo.DWG:

WMS 1.3 defines a "Map CS" (CRS:1) in Annex B.2, with origin in upper left. I believe we added that based on comments during the ISO 19128 review phase stating that we had not adequately defined how the geographic projection corresponds to the map view.

[Alex Dumitru]

@Dimitar: That is just a convention on the WMS service, they like to restrict arbitrary things for no reason at all

I looked over the topic yesterday and I believe this is the *only* correct solution. In pseudo-python below:

  def get_origin:
     if geotransform[1] > 0: # in our concepts geotransfor[3] is x-axis offset vector
        origin_on_x_axis = geotransform[0] # 0 position on x axis
     else:
        origin_on_y_axis = geotransform[0] + raster_pixels_on_x * geotransform[3] # max position on x axis

     if geotransform[5] > 0: # in our concepts geotransform[3] is y-axis offset vector
        origin_on_x_axis = geotransform[3] # 0 position on y axis
     else:
        origin_on_y_axis = geotransform[3] + raster_pixels_on_y * geotransform[6] # max position on y axis

     return  (origin_on_x_axis, origin_on_y_axis)

  def get_extents:
     origin_x, origin_y = get_origin()
     extent_x = origin_x + raster_pixels_on_x * geotransform[1]
     extent_y = origin_y + raster_pixels_on_y * geotransofrm[5]
     return ([origin_x, extent_x], [origin_y, extent_y])     

[Bang Pham Huu]

@Alex: one question (without verifying your code), the result is top-left or lower-left with RGB when using wcst-import? same question with import georeferenced image when using wcst-import as you said this is the 'only' correct solution.

It looks like you have a small mistake in here I think you mean geotransform[5]

geotransform[6] # max position on y axis

What I want to see is import image by wcst-import or rasql in test services will have the same origin, so if your 'only solution' is right, then we have 2 cases in here:

+ wcst_import and rasql import same the lower left origin (it is opposite as Meteo.DWG - we still need to reconfirm).

+ wcst_import and rasql is different origin (wcst_import is upper-left and rasql is lower-left), the question is still here but as you said you are right then rasql need to change (and that is what I expect).

Please note that I do not claim anything with you in here, I just want to describe my idea is:

+ Both cases need the same results.
+ When I query from lowest point suppose X (0: 400), Y (0: 400) then I want to see:

- With kind of image that upper-left (Y1) not > upper-right (Y2), 
it should be bottom up.

- With kind of image that upper-left (Y1 = 0) < upper-right (Y2),
 it should be top down.

[Alex Dumitru]

Replying to bphamhuu:

@Alex: one question (without verifying your code), the result is top-left or lower-left with RGB when using wcst-import? same question with import georeferenced image when using wcst-import as you said this is the 'only' correct solution.

As the geotransform of the anthur.tif image is (0.0, 1.0, 0.0, 0.0, 0.0, 1.0), it means that the image will have the origin in the lower bottom corner.

It looks like you have a small mistake in here I think you mean geotransform[5]

geotransform[6] # max position on y axis

Yes, you are right, it should be 5 in there.

What I want to see is import image by wcst-import or rasql in test services will have the same origin, so if your 'only solution' is right, then we have 2 cases in here:

+ wcst_import and rasql import same the lower left origin (it is opposite as Meteo.DWG - we still need to reconfirm).

This will be the case. The Meteo.DWG decision does not apply to coverages, it applies to layers in WMS. Furthermore this is just an internal decision of them, it's not even standardized in WMS. The coverage ecosystem, fortunately, is flexible enough to avoid the need of deciding such things.

+ wcst_import and rasql is different origin (wcst_import is upper-left and rasql is lower-left), the question is still here but as you said you are right then rasql need to change (and that is what I expect).

Please note that I do not claim anything with you in here, I just want to describe my idea is:

+ Both cases need the same results.
+ When I query from lowest point suppose X (0: 400), Y (0: 400) then I want to see:

- With kind of image that upper-left (Y1) not > upper-right (Y2), 
it should be bottom up.

The answer to your questions comes from the DescribeCoverage request, you cannot determine that only by looking at the extent. In there the origin is specified which together with the offset vector would tell you what will happen when you subset using those coordinates. Let's take some examples to make it easier to understand:

Coverage 1: Origin=(1000, 1000), OffsetVector=([2, 0], [0, -1]), Raster=(0:400, 0:300) ⇒ Extent=(1000:1800, 700:1000) - this coverage has origin upper left of the extent
geo_subset([1000:1200, 700:800], "Coverage 1") ⇒ raster_subset([0:100, 200:300])

Coverage 2: Origin=(1000, 1000), OffsetVector=([2, 0], [0, 1]), Raster=(0:400, 0:300) ⇒ Extent=(1000:1800, 1000:1300) - this coverage has origin bottom left of the extent
geo_subset([1000:1200, 1000:1100], "Coverage 2") ⇒ raster_subset([0:100, 0:100])

Coverage 3: Origin=(1000, 1000), OffsetVector=([-1, 0], [0, 3]), Raster=(0:400, 0:300) ⇒ Extent=(600:1000, 1000:1900) - this coverage has origin bottom right of the extent
geo_subset([600:700, 1000:1100], "Coverage 3") ⇒ raster_subset([300:400, 0:100])

You might even have a coverage referenced by arrays (as opposed to offset vectors): Origin=(0, 5], RefArray=[(-1,7), (0,5), (0,4), (1,9)] Raster(0:1, 0:1), Extent=(-1:1, 5:9) - here as you can see the origin is not in any corner.

• With kind of image that upper-left (Y1 = 0) < upper-right (Y2), it should be top down.

}}}

As you can see, the coverage ecosystem is a lot more flexible than the usual GIS tools (which treat date as glorified pictures), so you have to be a bit careful before considering what the possible use-cases are to a problem. Hope it clears things up.

[Bang Pham Huu]

@Alex, thanks for example, it is more interesting to discuss. One question: what happen when the 'origin' in Coverage 3 bottom right of the extent, sorry for my bad because I could understand:

my query: select x(0: 400), y(0: 100) from a (a is (0: 400, 0: 400) )

+ Coverage 1: origin is upper left then I get a piece of top-down.
+ Coverage 2: origin is lower left then I get a piece of bottom-up.

What is about coverage 3? I hope I still could get a piece of image or I will query outside of image based on your origin?

Second question: what is your solution for this ticket with the new name "Recognize if an image is geo-referenced or not with GDAL". It is best when you can have a conclude that with this kind of origin then you can decide which image is non-georeference and georeference (mine is I could not do this, so as you can understand 'coverage ecosystem' in flexible then I would like to know how to solve this).

Please note that, I asked because I still don't see the big picture in here, your explanation is good enough but I need a conclusion from my questions. Hope you could help.

[Peter Baumann]

I maintain that from an image file that does not carry an identifiable CRS we cannot deduce a georeference. Guesswork on coordinate constellations does not work. All the coverage cases above can be deduced from a CRS _if_ there is a CRS. top or bottom left does not say anything.
So, no CRS ⇒ no georeference ⇒ make it a GridCoverage using IndexCRS.

[Alex Dumitru]

Replying to bphamhuu:

@Alex, thanks for example, it is more interesting to discuss. One question: what happen when the 'origin' in Coverage 3 bottom right of the extent, sorry for my bad because I could understand:

my query: select x(0: 400), y(0: 100) from a (a is (0: 400, 0: 400) )

As I said above, you cannot subset a coverage without knowing the origin, offset vector (or the position array as in Coverage 4) and the raster extent. The geo extent does not matter. What are these 3 properties of a? If given I can answer the query, if not it's not possible. Also, you are mixing rasql and wcps and it makes it hard to understand if the query refers to pixels or geopixels.

+ Coverage 1: origin is upper left then I get a piece of top-down.
+ Coverage 2: origin is lower left then I get a piece of bottom-up.

You might be mixing up the geo-axes with rasql domains. They do not correspond 1 to 1 in neither position or values. x(0:400) is not defined in any of the coverages in my previous post.

What is about coverage 3? I hope I still could get a piece of image or I will query outside of image based on your origin?

Second question: what is your solution for this ticket with the new name "Recognize if an image is geo-referenced or not with GDAL". It is best when you can have a conclude that with this kind of origin then you can decide which image is non-georeference and georeference (mine is I could not do this, so as you can understand 'coverage ecosystem' in flexible then I would like to know how to solve this).

comment:34 explains the algorithm that fixes the behavior.

Please note that, I asked because I still don't see the big picture in here, your explanation is good enough but I need a conclusion from my questions. Hope you could help.

[Bang Pham Huu]

Hi Alex, I would like to confirm your solution, after that I can learn better from what you have proved. Please let me know, what is the value of:

+ raster_pixels_on_x and + raster_pixels_on_y

(I suppose it is the width and height of a coverage in pixels ?).

[Bang Pham Huu]

After discussion yesterday, I can see that Vlad and Alex has their points in using the direction of offset vector in choosing origin of coverage, so Alex, please make a patch with wcst-import tool when you can apply your solution and after that we can see the same result when query by WCPS (now WCS can return PNG also (wait for patch being applied - ticket 876 ) and rasql.

[Bang Pham Huu]

As my ticket has been closed in here ​http://www.rasdaman.org/ticket/1018 due to duplicate with this ticket, it is around 1 month and the problem still stays in there when using wcst_import with non-georeferenced images like RBG, mr,…So I would like to ask Alex solve the ticket 1018 first (in fact, I've done like a quick fix for this but this is not elegant from your view) and with another special cases you can make a complex solution when you have done completely (I suppose you need more time to solve). Thanks.

Service Call: http://localhost:8080/rasdaman/ows?service=WCS&version=2.0.1&request=UpdateCoverage&coverageId=rgbx3&subset1=i(0.0,400.0)&inputCoverageRef=file:///tmp/b72c0cff_3b5c_4cbe_ac14_0c4a66700b3c.gml&subset2=j(344.0,0.0)
Error Code: RuntimeError
Error Text: <?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<ows:ExceptionReport version="2.0.0"
    xsd:schemaLocation="http://www.opengis.net/ows/2.0 http://schemas.opengis.net/ows/2.0/owsExceptionReport.xsd"
    xmlns:ows="http://www.opengis.net/ows/2.0" xmlns:xsd="http://www.w3.org/2001/XMLSchema-instance" xmlns:xlink="http://www.w3.org/1999/xlink">
    <ows:Exception exceptionCode="RuntimeError">
        <ows:ExceptionText>Runtime error while processing request</ows:ExceptionText>
    </ows:Exception>
    <ows:Exception exceptionCode="RuntimeError">
        <ows:ExceptionText>petascope.wcps2.error.managed.processing.UnorderedSubsetException : Invalid subsetting coordinates: 344.0:0.0 for axis j.</ows:ExceptionText>
    </ows:Exception>

</ows:ExceptionReport>

[Dimitar Misev]

What's the status on this ticket?

[Alex Dumitru]

Still open, unfortunately.
@Vlad: Can you please check the mapping from the geo-origin to the grid origin for all Coverage types and report back on this ticket? It seems to me that there is a discrepancy between Grid and RectifiedGrid coverage types regarding this mapping.

[Bang Pham Huu]

When fixing ticket #1244, I think maybe it can relate to this ticket, then I gave a try and actually it can be the answer. The problem is in WCPS, it does not convert to "CRS:1" in case of grid axis. I've fixed it in WCPS 1.5 and now they has the same origin (they are the coverages which were imported in "WCPS system test": rgb and "WCST_Import" rgbx12).

in WCPS 1.0, if explicitly set the "CRS:1" to axis then the result is same:

'for c in ('''rgb''') return encode(c[i:"CRS:1"(0:400), j:"CRS:1"(0:100)], "png")'

'for c in ('''rgbx12''') return encode(c[i:"CRS:1"(0:400), j:"CRS:1"(0:100)], "png")'

Below is the result which is in WCPS 1.5

https://drive.google.com/file/d/0B_VL2-SVlKPHT1VxVkVhUGlpelU/view?usp=sharing

@Vlad: I think (If I'm not wrong), this ticket can be closed as it is problem in WCPS without setting correct CRS Index for grid axis.

[Vlad Merticariu]

Ok Bang, good find, thanks!