Re: [gdal-dev] gdalwarp and NOAA L1b

[Andrew Brooks]

On Thu, 14 May 2009 20:21:11 +0100, yehya imam  wrote:
>
> I've downloaded AVHRR L1b images from NOAA CLASS. I'm trying to rectify the 
> images via gdalwarp -tps. The output is obviously erroneous as per gdalinfo 

The error message "There is a problem to invert the interpolation matrix" 
explains why gdalinfo gave erroneous information.

I too have had great difficulties with gdalwarp and L1B.

The way the driver reports the geolocation data is not optimal and the methods 
used to process the data are also not optimal so it's no surprise that the 
result is often broken.

I submitted a patch to fix one of the issues but it was never applied; you can 
download it from here: http://trac.osgeo.org/gdal/ticket/2501

I also submitted some code to perform proper interpolation but nobody has made 
use of it.

I started a dialogue with the driver author who said he was going to rewrite it 
but then he went silent and I've not heard from him since.

One method I've used in the past to get some results is as follows:

1) Extract ALL of the geolocation data myself (51 points per scan line) since 
the driver normally subsamples
2) Interpolate this myself (Lagrangian method) to get 2048 points per scan line
3) Create two VRT files for the latitude and longitude data
4) Create a VRT file describing the image data in the file which references the 
two geolocation VRT files
5) Use the VRT file with my patched version of gdalwarp

This was successful for me in the past but when I tried it again recently it 
failed.  I can give you more details if you want though.

Andrew

___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] gdalwarp and NOAA L1b

[Markus Neteler]

On Fri, May 15, 2009 at 1:46 PM, Andrew Brooks  wrote:
> On Thu, 14 May 2009 20:21:11 +0100, yehya imam  wrote:
>>
>> I've downloaded AVHRR L1b images from NOAA CLASS. I'm trying to rectify the 
>> images via gdalwarp -tps. The output is obviously erroneous as per gdalinfo
>
> The error message "There is a problem to invert the interpolation matrix" 
> explains why gdalinfo gave erroneous information.
>
> I too have had great difficulties with gdalwarp and L1B.

Me, too. I gave up so far but your answer raised my interest. I still
want to process large amounts of NOAA L1b.

> The way the driver reports the geolocation data is not optimal and the 
> methods used to process the data are also not optimal so it's no surprise 
> that the result is often broken.
>
> I submitted a patch to fix one of the issues but it was never applied; you 
> can download it from here: http://trac.osgeo.org/gdal/ticket/2501

(sidenote: I see you aren't in CC in that ticket)

> I also submitted some code to perform proper interpolation but nobody has 
> made use of it.

Was that submitted as ticket, too? I would like to try it.

> I started a dialogue with the driver author who said he was going to rewrite 
> it but then he went silent and I've not heard from him since.
>
> One method I've used in the past to get some results is as follows:
>
> 1) Extract ALL of the geolocation data myself (51 points per scan line) since 
> the driver normally subsamples
> 2) Interpolate this myself (Lagrangian method) to get 2048 points per scan 
> line
> 3) Create two VRT files for the latitude and longitude data
> 4) Create a VRT file describing the image data in the file which references 
> the two geolocation VRT files
> 5) Use the VRT file with my patched version of gdalwarp
>
> This was successful for me in the past but when I tried it again recently it 
> failed.  I can give you more details if you want though.

Perhaps it is worth to document that in the trac Wiki?

thanks
Markus
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] about SetRasterColorInterpretation

[Frank Warmerdam]

jor sion wrote:

Hi,
 
I used the the SetRasterColorInterpretation() method to modified the 
ColorInterpretation of a band, but it failed.
 
so SetRasterColorInterpretation() doesn't work?


If I want to modify the ColorInterpretation of a band, what should I do?


JoSn,

This method is not operational for very many drivers.

I see that if the PAM mechanism is enabled all PAM enabled drivers will
actually store the color interpretation in the .aux.xml file.  What version
of GDAL are you using?  What file format?

Best regards,
--
---+--
I set the clouds in motion - turn up   | Frank Warmerdam, warmer...@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush| Geospatial Programmer for Rent

___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] (no subject)

[Frank Warmerdam]

Sjur Kolberg wrote:
 
Hello, GDAL list;
 
Can anyone suggest what prevents the following GDALSuggestedWarpOutput 
from calculating an approximate geotransform? The resulting suggested 
GeoTransform is just a copy of the original GeoTransform.
 
Version is GDAL 1.6.0beta1 under Windows XP, I use Visual Studio 2005, 
this is a MFC application. poDataset is a 9900 band GRIB 1 data set 
from ECMWF. It appears to be read in OK, I can extract 
plausible metadata both from the dataset and from each band. 
 
   SourceWKT is (CString)

GEOGCS["Coordinate System imported from GRIB file",DATUM["unknown",
SPHEROID["Sphere",6367470,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]
   SourceWKT is from poDataset->GetProjectionRef()
 
   TargetWKT is (CString)

PROJCS["UTM 32 North (WGS84)",GEOGCS["WGS 84",DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],
TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",
0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9108"]],AXIS["Lat",NORTH],AXIS["Long",EAST],
AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],
PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",50],
PARAMETER["false_northing",0],UNIT["Meter",1]]
   TargetWKT is from 
ogrspatref.SetProjCS(msg);// CString msg is just a name

ogrspatref.SetWellKnownGeogCS( "WGS84" );
if (ogrspatref.SetUTM( zone, i ) != OGRERR_NONE)  // zone==32, i==1 for 
northern hem.
 
 
And the code goes:



void *hTransformArg = GDALCreateGenImgProjTransformer( 
poDataset, sourceWKT, NULL, targetWKT, FALSE, 0, 1 );


CPLAssert( hTransformArg != NULL );

if (GDALSuggestedWarpOutput(pSrcDS, GDALGenImgProjTransform, 
hTransformArg, SugGeoTransform, &nPixels, &nLines )  != CE_None )


Sjur,

I tried this (test.wkt is the UTM WKT) and it seems to work plausibly:

testepsg -t 'GEOGCS["Coordinate System imported from GRIB 
file",DATUM["unknown",SPHEROID["Sphere",6367470,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]' 
test.wkt 4 65

OGRCT: Source: +proj=longlat +a=6367470 +b=6367470 +no_defs
OGRCT: Target: +proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
(4.00,65.00,0.00) -> (264409.202311,7217779.091236,0.00)

So, the reprojection operation can apparently work.  The behavior you
are seeing appears to be congruent with reprojection doing nothing.
Is there any possibility that the PROJ library is not built into your
GDAL?

I see no other reason things should not work as you have constructed
them.  All I can see is to either debug through and see what is going
wrong, or to submit a ticket with info on how to reproduce the problem.

Best regards,
--
---+--
I set the clouds in motion - turn up   | Frank Warmerdam, warmer...@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush| Geospatial Programmer for Rent

___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] GDALWarp GCP

[Matt Hanson]

Hello,   I sent the email to the list in April but it's gotten no
response.I've since scoured the 'net and the GDAL archives, tried
countless number of permutations and every time end up with a blank
image.   I'm obviously missing something, yet the archives and
documentation don't provide any clear explanation of how this may be
accomplished.Any pointers would be greatly appreciated!
Thanks!
Matt

Wed, 22 Apr 2009 08:28:58 -0700
Hello, 

I'm having trouble getting an image warped  properly using GCPs.   I'm
wondering if I'm missing a piece in my code.

Here  I am warping from an image with a presumably unknown projection to
a georeferenced reference image.   The images have been registered
either manually or automatically and the GCPs come from that -
pixel/line in the unknown image corresponding to georeferenced
coordinates in the reference image.


GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();

psWarpOptions->pTransformerArg = GDALCreateGCPTransformer(gdalgcpssize,
gdalgcps, order, 0);
psWarpOptions->pfnTransformer = GDALGCPTransform;
psWarpOptions->eResampleAlg = GRA_NearestNeighbour;

psWarpOptions->hSrcDS = // Set to source image
psWarpOptions->nBandCount = 0; 
psWarpOptions->pfnProgress = GDALTermProgress;   

double adfDstGeoTransform[6];
int nPixels=0, nLines=0;
GDALSuggestedWarpOutput( psWarpOptions->hSrcDS, GDALGCPTransform, 
psWarpOptions->pTransformerArg, adfDstGeoTransform, &nPixels, &nLines


// Create new output image of size nPixels x nLines
// Copy Projection info from georeferenced image
psWarpOptions->hDstDS = // Set to new output image

// Initialize and execute the warp operation
GDALWarpOperation oOperation;
oOperation.Initialize( psWarpOptions );
oOperation.ChunkAndWarpImage( 0, 0, nPixels, nLines );


Are these essentially the correct steps to do this?GDAL runs without
any errors and generates an output image which has the correct
projection info and reasonable looking geographic corner coordinates.
However the resulting image is all black.

Thanks in advance for any clues anyone could provide!

matt


___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] GDALWarp GCP

[Even Rouault]

Matt,

The symptoms you get are compatible with a common mistake : forgetting to 
close the output dataset with GDALClose().

Otherwise, if it's a more serious issue, I'd advise you studying carefully the 
source code of gdalwarp : 
http://trac.osgeo.org/gdal/browser/trunk/gdal/apps/gdalwarp.cpp

(I've tested that if I modify gdalwarp.cpp to remove the closing of the output 
dataset, I just get your symptoms)

Best regards,

Even
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] [ANN] GSDView v. 0.5.9

[Antonio Valentino]

Hi list,
I'm pleased to announce the availability of GSDView v. 0.5.9.

Geo-Spatial Data Viewer (GSDView) is a lightweight viewer for
geo-spatial data and products.
It is written in python and Qt4 and uses the GDAL library for data
access.

GSDView is modular and has a simple plug-in architecture.

You can find it at

http://sourceforge.net/projects/gsdview


-- 
Antonio Valentino
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] [ANN] GSDView v. 0.5.9

[Antonio Valentino]

Il giorno Fri, 15 May 2009 14:25:25 -0400
Greg Coats  ha scritto:

> Does GSDView support displaying GeoTIFF and GeoJPEG2000 images?
> Greg

GSDView should be able to open all file formats supported by the local
installation of GDAL.

Currently it is not able to handle RGB images and raster bands with
sighed data types. 

In case of raster bands with complex data type the module is extracted.

Best regards 

-- 
Antonio Valentino
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] [ANN] GSDView v. 0.5.9

[Antonio Valentino]

Il giorno Fri, 15 May 2009 14:49:52 -0400
Greg Coats  ha scritto:

> All of my images are color images with 8 bits for Red, 8 bits for  
> Green, 8 bits for Blue.
> You seem to be saying GSDView can display all GDAL supported image  
> formats, but that GSDView can only display grey scale images.

You can *open* a dataset containing an RGB image but can't display it
in RGB mode. You can only display each single 8bit raster band in
gray scale mode.

Visualization of colour images will be supported in future versions of
the SW.


-- 
Antonio Valentino
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Motion: Promote GDAL 1.6.1 RC2 to final

[Howard Butler]

Folks,

I declare this motion passed.  I will rename the RC2 file and post a  
message to -announce shortly.


Thanks,

Howard
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] GDAL 1.6.1 Released

[Howard Butler]

The GDAL team is pleased to announce the release of GDAL 1.6.1 final.   
GDAL 1.6.1 is the first maintenance release of the 1.6 branch, and it  
includes fixes for nearly 120 bugs.  See the release notes at  for more details.


- http://download.osgeo.org/gdal/gdal-1.6.1.tar.gz
- http://download.osgeo.org/gdal/gdal161.zip

Howard
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

RE: [gdal-dev] (no subject) > GDALSuggestedWarpOutput() only returns the unaltered source geotransform

[Sjur Kolberg]

Yes, thanks, Frank. That was it.

Or even easier, I simply copied a proj.dll from another location into the 
directory where my program (and gdal.dll) resides.

Is it so that gdal looks for proj.dll without being linked with proj.lib, and 
without complaining if it's not found?

I'll build a gdal with debug info, allowing me to track.

Best regards and thanks,

Sjur :-)


 

-Original Message-
From: Frank Warmerdam [mailto:warmer...@pobox.com] 
Sent: Friday, May 15, 2009 4:34 PM
To: Sjur Kolberg
Cc: 'gdal-dev@lists.osgeo.org'
Subject: Re: [gdal-dev] (no subject)


Sjur,

I tried this (test.wkt is the UTM WKT) and it seems to work plausibly:

testepsg -t 'GEOGCS["Coordinate System imported from GRIB 
file",DATUM["unknown",SPHEROID["Sphere",6367470,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'
 
test.wkt 4 65
OGRCT: Source: +proj=longlat +a=6367470 +b=6367470 +no_defs
OGRCT: Target: +proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
(4.00,65.00,0.00) -> (264409.202311,7217779.091236,0.00)

So, the reprojection operation can apparently work.  The behavior you are 
seeing appears to be congruent with reprojection doing nothing.
Is there any possibility that the PROJ library is not built into your GDAL?

I see no other reason things should not work as you have constructed them.  All 
I can see is to either debug through and see what is going wrong, or to submit 
a ticket with info on how to reproduce the problem.

Best regards,
-- 
---+
---+--
I set the clouds in motion - turn up   | Frank Warmerdam, warmer...@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush| Geospatial Programmer for Rent



Sjur Kolberg wrote:
>  
> Hello, GDAL list;
>  
> Can anyone suggest what prevents the following GDALSuggestedWarpOutput 
> from calculating an approximate geotransform? The resulting suggested 
> GeoTransform is just a copy of the original GeoTransform.
>  
> Version is GDAL 1.6.0beta1 under Windows XP, I use Visual Studio 2005, 
> this is a MFC application. poDataset is a 9900 band GRIB 1 data set 
> from ECMWF. It appears to be read in OK, I can extract plausible 
> metadata both from the dataset and from each band.
>  
>SourceWKT is (CString)
> GEOGCS["Coordinate System imported from GRIB file",DATUM["unknown", 
> SPHEROID["Sphere",6367470,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]
>SourceWKT is from poDataset->GetProjectionRef()
>  
>TargetWKT is (CString)
> PROJCS["UTM 32 North (WGS84)",GEOGCS["WGS 84",DATUM["WGS_1984", 
> SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],
> TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",
> 0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,
> AUTHORITY["EPSG","9108"]],AXIS["Lat",NORTH],AXIS["Long",EAST],
> AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],
> PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],
> PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",50],
> PARAMETER["false_northing",0],UNIT["Meter",1]]
>TargetWKT is from 
> ogrspatref.SetProjCS(msg);// CString msg is just a 
> name
> ogrspatref.SetWellKnownGeogCS( "WGS84" ); if (ogrspatref.SetUTM( zone, 
> i ) != OGRERR_NONE)  // zone==32, i==1 for northern hem.
>  
>  
> And the code goes:
> 
> 
> void *hTransformArg = GDALCreateGenImgProjTransformer( poDataset, 
> sourceWKT, NULL, targetWKT, FALSE, 0, 1 );
> 
> CPLAssert( hTransformArg != NULL );
> 
> if (GDALSuggestedWarpOutput(pSrcDS, GDALGenImgProjTransform, 
> hTransformArg, SugGeoTransform, &nPixels, &nLines )  != CE_None )
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev