[gdal-dev] Geolocation arrays - location interpretation

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

[gdal-dev] Geolocation arrays - location interpretation

Agram, Piyush S (334D)

Hi,

     I’m trying to track down systematic pixel shift effects observed with warping data using geolocation arrays:

https://trac.osgeo.org/gdal/ticket/6959

 

Essentially, the same data (image and location) provided as geolocation arrays is warped to a different output depending on how the arrays are laid out – WESN, WENS, EWSN, EWNS.

Gdalwarp produces least shift in warped data when the input arrays themselves are oriented WENS.

 

Is there a definitive definition for interpretation of lat,lon,value when geolocation arrays – i.e, are the lat/lon representative of pixel center (which I think should be the intent)?

Does GDAL by default assume that the geolocation value corresponds to the top/left of the pixel? Could this confusion result in the observed systematic shift?

 

Alternately, could the input arrays be reoriented to WENS using VRTs instead of rewriting large raster files?

 

Piyush


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Geolocation arrays - location interpretation

Even Rouault-2

On jeudi 28 septembre 2017 02:14:53 CEST Agram, Piyush S (334D) wrote:

> Hi,

> I’m trying to track down systematic pixel shift effects observed with

> warping data using geolocation arrays:

> https://trac.osgeo.org/gdal/ticket/6959

>

> Essentially, the same data (image and location) provided as geolocation

> arrays is warped to a different output depending on how the arrays are laid

> out – WESN, WENS, EWSN, EWNS.

Gdalwarp produces least shift in warped data

> when the input arrays themselves are oriented WENS.

> Is there a definitive definition for interpretation of lat,lon,value when

> geolocation arrays – i.e, are the lat/lon representative of pixel center

> (which I think should be the intent)?

 

That's a good question. I would assume that center of pixel would be what is expected in the gdalgeoloc code. I don't see anything explicitly said about that in

https://trac.osgeo.org/gdal/wiki/rfc4_geolocate

 

During your examination of the code, could you find what was expected ?

 

I assume you are using geolocation transformer with one of the drivers that support reporting geolocation arrays ? One possibility would be an inconsistency in the convention this driver (or the datasource itself) reports the lon,lat values with what the geolocation transformer expects.

 

 

Does GDAL by default assume that the

> geolocation value corresponds to the top/left of the pixel? Could this

> confusion result in the observed systematic shift?

> Alternately, could the input arrays be reoriented to WENS using VRTs instead

> of rewriting large raster files?

 

X_DATASET and Y_DATASET can point to any valid GDAL dataset, so a VRT should be possible (and you'll need a VRT to modify the values of X_DATASET and Y_DATASET ;-))

 

Even

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Geolocation arrays - location interpretation

Even Rouault-2
In reply to this post by Agram, Piyush S (334D)

On jeudi 28 septembre 2017 02:14:53 CEST Agram, Piyush S (334D) wrote:

> Hi,

> I’m trying to track down systematic pixel shift effects observed with

> warping data using geolocation arrays:

> https://trac.osgeo.org/gdal/ticket/6959

>

> Essentially, the same data (image and location) provided as geolocation

> arrays is warped to a different output depending on how the arrays are laid

> out – WESN, WENS, EWSN, EWNS.

Gdalwarp produces least shift in warped data

> when the input arrays themselves are oriented WENS.

> Is there a definitive definition for interpretation of lat,lon,value when

> geolocation arrays – i.e, are the lat/lon representative of pixel center

> (which I think should be the intent)?

 

That's a good question. I would assume that center of pixel would be what is expected in the gdalgeoloc code. I don't see anything explicitly said about that in

https://trac.osgeo.org/gdal/wiki/rfc4_geolocate

 

During your examination of the code, could you find what was expected ?

 

I assume you are using geolocation transformer with one of the drivers that support reporting geolocation arrays ? One possibility would be an inconsistency in the convention this driver (or the datasource itself) reports the lon,lat values with what the geolocation transformer expects.

 

 

Does GDAL by default assume that the

> geolocation value corresponds to the top/left of the pixel? Could this

> confusion result in the observed systematic shift?

> Alternately, could the input arrays be reoriented to WENS using VRTs instead

> of rewriting large raster files?

 

X_DATASET and Y_DATASET can point to any valid GDAL dataset, so a VRT should be possible (and you'll need a VRT to modify the values of X_DATASET and Y_DATASET ;-))

 

Even

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Geolocation arrays - location interpretation

Agram, Piyush S (334D)

Hi Even,

     I think I have figured out how the GeoLocTransformer works. Dumped out the backmap and might have some changes.

 

1)       Rounding up to the nearest integer is definitely the issue here.

2)       I’m going to try this first – maybe a weighted approach will be better. For now, if a pixel falls in the middle of square – I’m going to add it to all of the 4 vertices. Then average all the data that went into a pixel of the backmap.

3)       For the regular grid tests, if this shows improvement (hopefully) – I can share the implementation and we can possibly discuss about adding a weighting scheme. I don’t want to increase the runtime too much by adding complexity.

 

Piyush

 

From: Even Rouault <[hidden email]>
Date: Thursday, September 28, 2017 at 2:56 AM
To: "[hidden email]" <[hidden email]>
Cc: "Agram, Piyush S (334D)" <[hidden email]>
Subject: Re: [gdal-dev] Geolocation arrays - location interpretation

 

On jeudi 28 septembre 2017 02:14:53 CEST Agram, Piyush S (334D) wrote:

> Hi,

> I’m trying to track down systematic pixel shift effects observed with

> warping data using geolocation arrays:

> https://trac.osgeo.org/gdal/ticket/6959

>

> Essentially, the same data (image and location) provided as geolocation

> arrays is warped to a different output depending on how the arrays are laid

> out – WESN, WENS, EWSN, EWNS.

Gdalwarp produces least shift in warped data

> when the input arrays themselves are oriented WENS.

> Is there a definitive definition for interpretation of lat,lon,value when

> geolocation arrays – i.e, are the lat/lon representative of pixel center

> (which I think should be the intent)?

 

That's a good question. I would assume that center of pixel would be what is expected in the gdalgeoloc code. I don't see anything explicitly said about that in

https://trac.osgeo.org/gdal/wiki/rfc4_geolocate

 

During your examination of the code, could you find what was expected ?

 

I assume you are using geolocation transformer with one of the drivers that support reporting geolocation arrays ? One possibility would be an inconsistency in the convention this driver (or the datasource itself) reports the lon,lat values with what the geolocation transformer expects.

 

 

Does GDAL by default assume that the

> geolocation value corresponds to the top/left of the pixel? Could this

> confusion result in the observed systematic shift?

> Alternately, could the input arrays be reoriented to WENS using VRTs instead

> of rewriting large raster files?

 

X_DATASET and Y_DATASET can point to any valid GDAL dataset, so a VRT should be possible (and you'll need a VRT to modify the values of X_DATASET and Y_DATASET ;-))

 

Even

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Geolocation arrays - location interpretation

Even Rouault-2

On jeudi 28 septembre 2017 17:26:52 CEST Agram, Piyush S (334D) wrote:

> Hi Even,

> I think I have figured out how the GeoLocTransformer works. Dumped out

> the backmap and might have some changes.

>

> 1) Rounding up to the nearest integer is definitely the issue here.

>

> 2) I’m going to try this first – maybe a weighted approach will be

> better. For now, if a pixel falls in the middle of square – I’m going to

> add it to all of the 4 vertices. Then average all the data that went into a

> pixel of the backmap.

 

Perhaps bilinear resampling so that the weights are proportionnal ?

 

> 3) For the regular grid tests, if this shows improvement (hopefully) –

> I can share the implementation and we can possibly discuss about adding a

> weighting scheme. I don’t want to increase the runtime too much by adding

> complexity.

 

I think (sub-pixel) accuracy is more important than speed. Speed improvements can be investigated later if needed.

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Geolocation arrays - location interpretation

Agram, Piyush S (334D)

Hi Even,

    I implemented the bilinear weighting like we discussed and that took out the systematic pixel shift from 1 to 0.5 (sign depending on orientation).

You can see the changes here:

https://github.com/piyushrpt/gdal/commit/b199cbae9ea15e31c787c1480c342bb84debf774

 

I also tried playing with the upsampling factor (1.3 to 4.0) and the discrepancy numbers were consistent after 1.5. Probably, the minimum needs to be sqrt(2) for square pixels.

Using the weighting fixed the truncation / rounding error and the results seem consistent when I move the output extent around.

 

The systematic 0.5 pixel shift suggests that there is still an issue. Is my interpretation that the GeoLocTransformer is supposed to operate on pixel center coordinates correct?

i.e – if my geoloc arrays were regular grid and I feed in the coordinates to the transformer, they return pixel centers – 0.5,0.5 etc.

I’m going to try implementing some round trip checks to see if that shows something?  

 

Piyush

 

From: Even Rouault <[hidden email]>
Date: Thursday, September 28, 2017 at 12:25 PM
To: "Agram, Piyush S (334D)" <[hidden email]>
Cc: "[hidden email]" <[hidden email]>
Subject: Re: [gdal-dev] Geolocation arrays - location interpretation

 

On jeudi 28 septembre 2017 17:26:52 CEST Agram, Piyush S (334D) wrote:

> Hi Even,

> I think I have figured out how the GeoLocTransformer works. Dumped out

> the backmap and might have some changes.

>

> 1) Rounding up to the nearest integer is definitely the issue here.

>

> 2) I’m going to try this first – maybe a weighted approach will be

> better. For now, if a pixel falls in the middle of square – I’m going to

> add it to all of the 4 vertices. Then average all the data that went into a

> pixel of the backmap.

 

Perhaps bilinear resampling so that the weights are proportionnal ?

 

> 3) For the regular grid tests, if this shows improvement (hopefully) –

> I can share the implementation and we can possibly discuss about adding a

> weighting scheme. I don’t want to increase the runtime too much by adding

> complexity.

 

I think (sub-pixel) accuracy is more important than speed. Speed improvements can be investigated later if needed.

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Geolocation arrays - location interpretation

Even Rouault-2

On jeudi 28 septembre 2017 22:18:31 CEST Agram, Piyush S (334D) wrote:

> Hi Even,

> I implemented the bilinear weighting like we discussed and that took out

> the systematic pixel shift from 1 to 0.5 (sign depending on orientation).

> You can see the changes here:

> https://github.com/piyushrpt/gdal/commit/b199cbae9ea15e31c787c1480c342bb84de

> bf774

> I also tried playing with the upsampling factor (1.3 to 4.0) and the

> discrepancy numbers were consistent after 1.5. Probably, the minimum needs

> to be sqrt(2) for square pixels.

Using the weighting fixed the truncation

> / rounding error and the results seem consistent when I move the output

> extent around.

> The systematic 0.5 pixel shift suggests that there is still an issue. Is my

> interpretation that the GeoLocTransformer is supposed to operate on pixel

> center coordinates correct?

i.e – if my geoloc arrays were regular grid

> and I feed in the coordinates to the transformer, they return pixel centers

> – 0.5,0.5 etc. I’m going to try implementing some round trip checks to see

> if that shows something?

 

Yes I agree there must be a 0.5 pixel shift issue when looking at the forward path in GDALGeoLocTransform(), ie the one taken by if( !bDstToSrc )

 

So we have currently

 

const double dfGeoLocPixel =

(padfX[i] - psTransform->dfPIXEL_OFFSET)

/ psTransform->dfPIXEL_STEP;

const double dfGeoLocLine =

(padfY[i] - psTransform->dfLINE_OFFSET)

/ psTransform->dfLINE_STEP;

 

 

For example take the case of offset = 0 and step = 1 to simplify

 

Imagine than padfX[0] = 0.5 and padfY[0] = 0.5 (so center of the top left pixel)

 

which will lead to iX = 0 and iY = 10

 

Now we are going to take a bilinear resampling at lines 788-800 (current trunk)

 

out_padfX[0] = 0.25 * (geoloc[0] + geoloc[1] + geoloc[width] + geoloc[width+1])

 

whereas it should be just out_padfX[0] = geoloc[0] . No resampling needed.

 

So perhaps something like the following ( substracting 0.5 ) would work ?

 

const double dfGeoLocPixel =

(padfX[i] - 0.5 - psTransform->dfPIXEL_OFFSET)

/ psTransform->dfPIXEL_STEP;

const double dfGeoLocLine =

(padfY[i] - 0.5 - psTransform->dfLINE_OFFSET)

/ psTransform->dfLINE_STEP;

 

and probably something equivalent in the other direction

 

Perhaps we should have some metadata item to indicate the convention and apply or not this 0.5 shift

Like GEOLOC_COORD_CONVENTION=PIXEL_CENTER (default if not specified) or GEOLOC_COORD_CONVENTION=PIXEL_TOP_LEFT

But I'd expect most (all ?) datasources to use a pixel_center convention.

 

Which format/driver are you using ?

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Geolocation arrays - location interpretation

Agram, Piyush S (334D)

Hi Even,

     Here is my first cut at fixing the half pixel shifts:

https://github.com/piyushrpt/gdal/commit/ead79550e7120767187a0b9c836b5d861223d997

 

With these changes, the shifts are down to sub 0.2 pixels. Still not sure if I got the math right. Feel like the errors should be smaller for the ideal simulated datasets that I’m using.

 

Piyush

 

From: Even Rouault <[hidden email]>
Date: Thursday, September 28, 2017 at 3:36 PM
To: "Agram, Piyush S (334D)" <[hidden email]>
Cc: "[hidden email]" <[hidden email]>
Subject: Re: [gdal-dev] Geolocation arrays - location interpretation

 

On jeudi 28 septembre 2017 22:18:31 CEST Agram, Piyush S (334D) wrote:

> Hi Even,

> I implemented the bilinear weighting like we discussed and that took out

> the systematic pixel shift from 1 to 0.5 (sign depending on orientation).

> You can see the changes here:

> https://github.com/piyushrpt/gdal/commit/b199cbae9ea15e31c787c1480c342bb84de

> bf774

> I also tried playing with the upsampling factor (1.3 to 4.0) and the

> discrepancy numbers were consistent after 1.5. Probably, the minimum needs

> to be sqrt(2) for square pixels.

Using the weighting fixed the truncation

> / rounding error and the results seem consistent when I move the output

> extent around.

> The systematic 0.5 pixel shift suggests that there is still an issue. Is my

> interpretation that the GeoLocTransformer is supposed to operate on pixel

> center coordinates correct?

i.e – if my geoloc arrays were regular grid

> and I feed in the coordinates to the transformer, they return pixel centers

> – 0.5,0.5 etc. I’m going to try implementing some round trip checks to see

> if that shows something?

 

Yes I agree there must be a 0.5 pixel shift issue when looking at the forward path in GDALGeoLocTransform(), ie the one taken by if( !bDstToSrc )

 

So we have currently

 

const double dfGeoLocPixel =

(padfX[i] - psTransform->dfPIXEL_OFFSET)

/ psTransform->dfPIXEL_STEP;

const double dfGeoLocLine =

(padfY[i] - psTransform->dfLINE_OFFSET)

/ psTransform->dfLINE_STEP;

 

 

For example take the case of offset = 0 and step = 1 to simplify

 

Imagine than padfX[0] = 0.5 and padfY[0] = 0.5 (so center of the top left pixel)

 

which will lead to iX = 0 and iY = 10

 

Now we are going to take a bilinear resampling at lines 788-800 (current trunk)

 

out_padfX[0] = 0.25 * (geoloc[0] + geoloc[1] + geoloc[width] + geoloc[width+1])

 

whereas it should be just out_padfX[0] = geoloc[0] . No resampling needed.

 

So perhaps something like the following ( substracting 0.5 ) would work ?

 

const double dfGeoLocPixel =

(padfX[i] - 0.5 - psTransform->dfPIXEL_OFFSET)

/ psTransform->dfPIXEL_STEP;

const double dfGeoLocLine =

(padfY[i] - 0.5 - psTransform->dfLINE_OFFSET)

/ psTransform->dfLINE_STEP;

 

and probably something equivalent in the other direction

 

Perhaps we should have some metadata item to indicate the convention and apply or not this 0.5 shift

Like GEOLOC_COORD_CONVENTION=PIXEL_CENTER (default if not specified) or GEOLOC_COORD_CONVENTION=PIXEL_TOP_LEFT

But I'd expect most (all ?) datasources to use a pixel_center convention.

 

Which format/driver are you using ?

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Geolocation arrays - location interpretation

Even Rouault-2

On vendredi 29 septembre 2017 03:35:55 CEST Agram, Piyush S (334D) wrote:

> Hi Even,

> Here is my first cut at fixing the half pixel shifts:

> https://github.com/piyushrpt/gdal/commit/ead79550e7120767187a0b9c836b5d86122

> 3d997

> With these changes, the shifts are down to sub 0.2 pixels. Still not sure if

> I got the math right. Feel like the errors should be smaller for the ideal

> simulated datasets that I’m using.

 

Yes I think with a synthetic geolocation array that would encode the equivalent of a geotransform matrix (padfX[y][x] = X0 + x_res * x + 0.5 and padfY[y][x] = Y0 + y_res * y - 0.5), you should ideally get a ~0 pixel error.

One potential candidate for error would be the "loop over the backmap trying to fill in holes with nearby values." step. Given the backmap is created at a higher resolution than the source, the interpolation done in it can contribute errors.

On the above synthetic geolocation array, I would probably try with a backmap size equal to the geolocation dataset size (which should eliminate any hole filling step), and then tune things until I get a 0 error.

Perhaps the hole filling step would need a smarter/smoother interpolation logic

 

> Piyush

>

> From: Even Rouault <[hidden email]>

> Date: Thursday, September 28, 2017 at 3:36 PM

> To: "Agram, Piyush S (334D)" <[hidden email]>

> Cc: "[hidden email]" <[hidden email]>

> Subject: Re: [gdal-dev] Geolocation arrays - location interpretation

>

>

> On jeudi 28 septembre 2017 22:18:31 CEST Agram, Piyush S (334D) wrote:

>

>

> > Hi Even,

>

>

>

> > I implemented the bilinear weighting like we discussed and that took out

>

>

>

> > the systematic pixel shift from 1 to 0.5 (sign depending on orientation).

>

>

>

> > You can see the changes here:

>

>

>

> > https://github.com/piyushrpt/gdal/commit/b199cbae9ea15e31c787c1480c342bb84

> > de

>

>

> > bf774

>

>

>

> > I also tried playing with the upsampling factor (1.3 to 4.0) and the

>

>

>

> > discrepancy numbers were consistent after 1.5. Probably, the minimum

> > needs

>

>

>

> > to be sqrt(2) for square pixels.

>

>

> Using the weighting fixed the truncation

>

>

> > / rounding error and the results seem consistent when I move the output

>

>

>

> > extent around.

>

>

>

> > The systematic 0.5 pixel shift suggests that there is still an issue. Is

> > my

>

>

> > interpretation that the GeoLocTransformer is supposed to operate on pixel

>

>

>

> > center coordinates correct?

>

>

> i.e – if my geoloc arrays were regular grid

>

>

> > and I feed in the coordinates to the transformer, they return pixel

> > centers

>

>

> > – 0.5,0.5 etc. I’m going to try implementing some round trip checks to

> > see

>

>

>

> > if that shows something?

>

>

>

>

> Yes I agree there must be a 0.5 pixel shift issue when looking at the

> forward path in GDALGeoLocTransform(), ie the one taken by if( !bDstToSrc

> )

>

>

> So we have currently

>

>

>

> const double dfGeoLocPixel =

>

> (padfX[i] - psTransform->dfPIXEL_OFFSET)

>

> / psTransform->dfPIXEL_STEP;

>

> const double dfGeoLocLine =

>

> (padfY[i] - psTransform->dfLINE_OFFSET)

>

> / psTransform->dfLINE_STEP;

>

>

>

>

>

> For example take the case of offset = 0 and step = 1 to simplify

>

>

>

> Imagine than padfX[0] = 0.5 and padfY[0] = 0.5 (so center of the top left

> pixel)

>

>

> which will lead to iX = 0 and iY = 10

>

>

>

> Now we are going to take a bilinear resampling at lines 788-800 (current

> trunk)

>

>

> out_padfX[0] = 0.25 * (geoloc[0] + geoloc[1] + geoloc[width] +

> geoloc[width+1])

>

>

> whereas it should be just out_padfX[0] = geoloc[0] . No resampling needed.

>

>

>

> So perhaps something like the following ( substracting 0.5 ) would work ?

>

>

>

> const double dfGeoLocPixel =

>

> (padfX[i] - 0.5 - psTransform->dfPIXEL_OFFSET)

>

> / psTransform->dfPIXEL_STEP;

>

> const double dfGeoLocLine =

>

> (padfY[i] - 0.5 - psTransform->dfLINE_OFFSET)

>

> / psTransform->dfLINE_STEP;

>

>

>

> and probably something equivalent in the other direction

>

>

>

> Perhaps we should have some metadata item to indicate the convention and

> apply or not this 0.5 shift

> Like GEOLOC_COORD_CONVENTION=PIXEL_CENTER (default if not specified) or

> GEOLOC_COORD_CONVENTION=PIXEL_TOP_LEFT

> But I'd expect most (all ?) datasources to use a pixel_center convention.

>

>

>

> Which format/driver are you using ?

>

>

>

> --

>

> Spatialys - Geospatial professional services

>

> http://www.spatialys.com

 

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev