[postgis-devel] [gdal-dev] Question about gdal2wktraster

[Jorge Arévalo]

Hello,

On Tue, Jul 14, 2009 at 5:21 PM, Mateusz Loskot<mateusz at loskot.net> wrote:
> Jorge Arévalo wrote:
>>
>> Hello,
>>
>> On Tue, Jul 14, 2009 at 3:28 PM, Mateusz Loskot<mateusz at loskot.net> wrote:
>>>>>
>>>>> I'd suggest to test datasets based on EPSG 4326 first.
>>>>> After these work for you, then switch to projected ones.
>>>>
>>>> Do you know where can I find geotiff images based on EPSG:4326?
>>>
>>> Jorge,
>>>
>>> My favorite images for initial and basic testing of WKT Raster are here:
>>>
>>> http://download.osgeo.org/gdal/data/gtiff/
>>>
>>> First, I use utm.tif to test if 1 band raster works.
>>> This one is UTM (EPSG:26711) not WGS 84, but it's perfectly usable as
>>> such.
>>> Also, you can easily transform it to WGS84 with GDAL utils if you like.
>>
>> Actually, this is the same image I was using (cea.tif), but mine has
>> wrong srid... Thanks a lot for the info!
>
> Yes, I know. utm.tif and cea.tif have the same content
> but these datasets are different, not the same because of projection.

Ok, here are the Corner coordinates of the utm.tif file (gdalinfo output):

Corner Coordinates:
Upper Left  (  440720.000, 3751320.000) (117d38'28.21"W, 33d54'8.47"N)
Lower Left  (  440720.000, 3720600.000) (117d38'20.79"W, 33d37'31.04"N)
Upper Right (  471440.000, 3751320.000) (117d18'32.07"W, 33d54'13.08"N)
Lower Right (  471440.000, 3720600.000) (117d18'28.50"W, 33d37'35.61"N)
Center      (  456080.000, 3735960.000) (117d28'27.39"W, 33d45'52.46"N)

After loading the raster in PostGIS, the Corner Coordinates set by my
driver (gdalinfo ouptut):

Corner Coordinates:
Upper Left  (  440720.000, 3720600.000) (117d38'20.79"W, 33d37'31.04"N)
Lower Left  (  440720.000, 3689880.000) (117d38'13.48"W, 33d20'53.57"N)
Upper Right (  471440.000, 3720600.000) (117d18'28.50"W, 33d37'35.61"N)
Lower Right (  471440.000, 3689880.000) (117d18'24.98"W, 33d20'58.09"N)
Center      (  456080.000, 3705240.000) (117d28'21.94"W, 33d29'14.98"N)

As you can see, there are small differences in Y coordinate. I
populate the raster info in this way:

(...)
pszProjectionRef = (char *) GetProjectionRef(); // returns the correct
projection ref for srid 26711

poSR = new OGRSpatialReference(pszProjectionRef);

// pszWKTRaster has the extent of the raster, got by ST_AsText(extent)
from raster_columns
// poGeom will store the geometry
poOgrErr = OGRGeometryFactory::createFromWkt(&pszWKTRaster, poSR,
&poGeom);

poE = new OGREnvelope();
poGeom->getEnvelope(poE);

// raster coordinates
fUpperLeftX = poE->MinX;
fUpperLeftY = poE->MinY;

// fPixelSizeX and fPixelSizeY got from raster_columns
nRasterXSize = (int) fabs(rint((poE->MaxX - poE->MinX) / fPixelSizeX));
nRasterYSize = (int) fabs(rint((poE->MaxY - poE->MinY) / fPixelSizeY));
(...)

Then, I suppose that gdalinfo obtains corner coordinates from upper
left coordinates and raster size. Is there any obvious problem in this
code? Maybe the rint function?

Thanks in advance

Best regards
Jorge


>
> Best regards,
> --
> Mateusz Loskot, http://mateusz.loskot.net
> _______________________________________________
> postgis-devel mailing list
> postgis-devel at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-devel
>