[gdal-dev] reprojecting coastcolour (/meris) using python via GCPs

Tue Feb 18 05:15:46 PST 2014

On Tue, Feb 18, 2014 at 9:22 AM, Ivan Price <Ivan.Price at noveltis.fr> wrote:

>  Hi again Etienne,
>
>
>
> sorry for being vague.. when i say it 'works' in inverted commas, its
> because i'm not sure to what level it is expected to work. I can read the
> list of bands but there is no georeferencing for example., gdal doesn't see
> the gcps (=values in lat/lon variables).
>
>
>
> so previously i was writing out a (subsetted) gtiff with gcps in python
> then warping using gdalwarp -tps which was working well actually.. just a
> few more steps than i'd ideally like.
>
>
>
> following your advice regarding the vrt and after learning from nansat
> examples (thanks Anton!) i now write a vrt with gcps and use gdalwarp from
> the commandline with -tps on that.. which works fine.. although i need to
> use a gcp for every > 200 pixels or the calculation time is too high.
>
>
>
>
>
> if you tell me that the gcps should be recognised automatically by gdal
> then i'd be happy to file a bug and upload some sample data.
>
>
>
> for info i'm using gdal 1.9.2, built by me using flags:
>
> --prefix --with-python --pymoddir=xxx
>
> and thats all.
>

You might have better success using a more recent gdal version, the netcdf
driver has changed a bit since 1.9. I am not sure how it would handle gcps
though, I added support for GEOLOCATION arrays but that is a different
thing.

Seems to me using a vrt is the best option.

Etienne

>
>
> thanks for your assistance,
>
>
>
> -i
>
>
>
> save.nc is an untouched coastcolour dataset.
>
>
>
> (price at nemo:)$ gdalinfo NETCDF:"save.nc":reflec_1
>
> Driver: netCDF/Network Common Data Format
>
> Files: none associated
>
> Size is 4481, 9921
>
> Coordinate System is `'
>
> Metadata:
>
>   NC_GLOBAL#auto_grouping=reflec:norm_reflec:atm_tau
>
>   NC_GLOBAL#Conventions=CF-1.4
>
>   NC_GLOBAL#metadata_profile=beam
>
>   NC_GLOBAL#metadata_version=0.5
>
>   NC_GLOBAL#product_type=MER_FSG_CCL2R
>
>   NC_GLOBAL#start_date=08-FEB-2011 10:37:07.984543
>
>   NC_GLOBAL#stop_date=08-FEB-2011 10:44:24.434783
>
>   NC_GLOBAL#TileSize=16:4481
>
>   NC_GLOBAL#title=MERIS CoastColour L2R
>
>   reflec_1#bandwidth=9.9370003
>
>   reflec_1#coordinates=lat lon
>
>   reflec_1#long_name=Water leaving irradiance reflectance at 412.691 nm
>
>   reflec_1#solar_flux=1764.0343
>
>   reflec_1#spectral_band_index=0
>
>   reflec_1#units=sr^-1
>
>   reflec_1#valid_pixel_expression=!l2r_flags.INPUT_INVALID
>
>   reflec_1#wavelength=412.69101
>
> Corner Coordinates:
>
> Upper Left  (    0.0,    0.0)
>
> Lower Left  (    0.0, 9921.0)
>
> Upper Right ( 4481.0,    0.0)
>
> Lower Right ( 4481.0, 9921.0)
>
> Center      ( 2240.5, 4960.5)
>
> Band 1 Block=4481x1 Type=Float32, ColorInterp=Undefined
>
>   NoData Value=9.96920996838686905e+36
>
>   Metadata:
>
>     bandwidth=9.9370003
>
>     coordinates=lat lon
>
>     long_name=Water leaving irradiance reflectance at 412.691 nm
>
>     NETCDF_VARNAME=reflec_1
>
>     solar_flux=1764.0343
>
>     spectral_band_index=0
>
>     units=sr^-1
>
>     valid_pixel_expression=!l2r_flags.INPUT_INVALID
>
>     wavelength=412.69101
>
>
>
>
>
> (price at nemo:)$ gdalwarp -t_srs EPSG:4326  NETCDF:"save.nc":reflec_1
> foo.tif
>
> Processing input file NETCDF:save.nc:reflec_1.
>
> ERROR 1: Unable to compute a transformation between pixel/line
>
> and georeferenced coordinates for NETCDF:save.nc:reflec_1.
>
> There is no affine transformation and no GCPs.
>
>
>
>
>
>
>
>
>
> *De :* Etienne Tourigny [mailto:etourigny.dev at gmail.com]
> *Envoyé :* Monday, 17 February 2014 20:32
>
> *À :* Ivan Price
> *Cc :* gdal-dev at lists.osgeo.org
> *Objet :* Re: [gdal-dev] reprojecting coastcolour (/meris) using python
> via GCPs
>
>
>
>
>
>
>
> On Mon, Feb 17, 2014 at 2:03 PM, Ivan Price <Ivan.Price at noveltis.fr>
> wrote:
>
> Hi Etienne,
>
>
>
> firstly the problem reading the MERIS (actually coastcolour) data was due
> to the fact i was using gdal 1.6, after upgrading to 1.9 i can at least
> list the bands so i guess that 'works'.
>
>
>
> does it really work or sort of? If not, please file a bug report with an
> example data set attached or a link to it.
>
>
>
>
>
> regarding the use of the filesystem, i tried that, however it is slow to
> write out the whole dataset to disk, (i need to treat hundreds of images),
> so i was looking for a way to avoid this.. although it will be my fallback
> option if all else fails..
>
>
>
> have you tried using a .vrt file?
>
>
>
>
>
> either way i'll reply to the list with what we do in the end
>
>
>
> thanks for your response,
>
>
>
> -i
>
>
>
>
>
>
>
> *De :* Etienne Tourigny [mailto:etourigny.dev at gmail.com]
> *Envoyé :* Monday, 17 February 2014 14:09
> *À :* Ivan Price
> *Cc :* gdal-dev at lists.osgeo.org
>
>
> *Objet :* Re: [gdal-dev] reprojecting coastcolour (/meris) using python
> via GCPs
>
>
>
> As far as I know, the gdal warp api is not exposed to python.
>
>
>
> I have no idea on reading the MERIS data with the gdal netcdf driver -
> what I the problem?
>
>
>
> But you might be able to use this workaround: instead of creating a source
> dataset using the memory driver, create a file on disk (in gtiff format)
> with GCPS and write data to disk, and then use gdalwarp on that.
>
>
>
> Etienne
>
>
>
> On Mon, Feb 17, 2014 at 6:46 AM, Ivan Price <Ivan.Price at noveltis.fr>
> wrote:
>
>
>
> Hello,
>
>
>
> I am trying to reproject a window inside a coastcolour (=MERIS) image. As
> far as I can see GDAL cannot read the coastcolour data directly, so i am
> reading the coastcolour netcdf in python, building a source dataset using
> the memory driver, adding GCPS (1 for every 10th pixel) and writing the
> data to it, then reprojecting the source dataset to a destination dataset
> which is a spatial subset of the original in wgs84 lat/long.
>
>
>
> This works fine and is relatively fast, but the reprojection is not
> accurate, the results are out by about 6-10 pixels (in various directions).
> On reading the forums it seems if i was using gdalwarp i would be using
> -tps, however the ReprojectImage() function does not seem to offer this
> parameter ? And i don't have the option of using the commandline tool as
> even gdal 1.10 cannot recognise the coastcolour data.
>
>
>
> So i guess i have 2 questions.. has anyone had any success reading
> coastcolour data with the gdal command line tools, and secondly:
>
>
>
> how can i get ReprojectImage() to be more accurate, given i have a GCP for
> every pixel ?
>
>
>
> thanks and regards,
>
>
>
> -ivan
>
>
>
>
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140218/723a5542/attachment-0001.html>