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

Anton Korosov anton.korosov at nersc.no
Tue Feb 18 00:49:20 PST 2014


Hi Ivan,

if you are using Python, then options are usually given as a list of 
strings. I cannot find in the gdal docs at the moment (docs for python 
bindings are quite scarce), but you can check in Nansat:

https://github.com/nansencenter/nansat/blob/develop/vrt.py#L1641

So for your case it may be ['TPS=True'] or ['METHOD=GCP_TPS']

Regarding Nansat, we can add your mapper in the repository, or I can 
help a little with development of it (if you provide band description 
and an example file). Our mailing list is nansat-dev at googlegroups.com


Regards!
Anton

On 02/17/2014 05:51 PM, Ivan Price wrote:
> Hi Anton,
>
> thanks for responding,
>
> I have tried adding the tps option like this:
>
> gdal.ReprojectImage(input_ds, output_ds, input_ds.GetProjection(), output_ds.GetProjection(),
>                          gdal.GRA_NearestNeighbour, 0.5*1024*1024*1024, 0, None, {'tps': True})
>
> but it has no effect. despite the doco in the link you have there is a max of 9 arguments, (not 10), so i'm not sure its even being used for what it should be. (i'm using the python bindings from gdal 1.9)
>
> regarding nansat it looks very interesting..
>
> after an initial usage I see that i need to make a new mapper module for coastcolour data, presumably based on the MERIS L2 example.. as there are metadata fields that are missing in the coastcolour data that the envisat mapper is missing.
>
> it is interesting to see is that the reprojection is done in python, but using a VRT mechanism. maybe i can adapt this approach to get something working.
>
> thanks very much for the pointer
>
> -i
>
>
>
>
> -----Message d'origine-----
> De : gdal-dev-bounces at lists.osgeo.org [mailto:gdal-dev-bounces at lists.osgeo.org] De la part de Anton Korosov
> Envoyé : Monday, 17 February 2014 13:38
> À : gdal-dev at lists.osgeo.org
> Objet : Re: [gdal-dev] reprojecting coastcolour (/meris) using python via GCPs
>
> Hello Ivan,
>
> ReprojectImage() has a parameter psOptions()
> http://www.gdal.org/gdalwarper_8h.html#ad36462e8d5d34642df7f9ea1cfc2fec4
> It should accept any of the warping options, e.g. '-tps'
> http://www.gdal.org/gdalwarp.html
>
> I'm not 100% sure (since I'm using Python) but it may work.
>
>
> I'll take this chance also to promote Nansat, a scientist-friendly Python tool for working with satellite and model data:
> Repo: https://github.com/nansencenter/nansat
> Wiki: https://github.com/nansencenter/nansat/wiki
> and API-rference: http://nansencenter.github.io/nansat/
>
> It is a wrapper around GDAL, which adds scientific meaning to the opened images. Briefly: GDAL doesn't know much about e.g. band 22 in a MERIS image, Nansat does. It provides full information and allows simple usage, e.g.:
> n = Nansat(meris_image)
> n.reproject(dstDomain, tps=True)
> n.export('outFile.nc')
>
> It can open MERIS as well as tens of other formats.
>
>
> Best regards!
> Anton
>
> On 02/17/2014 10:46 AM, Ivan Price 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
>>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>



More information about the gdal-dev mailing list