[gdal-dev] transformation with nadgrids

2010-12-13 Thread Wolfgang Lahr
Hello list,

I am using GDAL/OGR to transform raster and vector data from German
Gauss-Krueger (DHDN) to ETRS/UTM like so:

gdalwarp/ogr2ogr -s_srs +proj=tmerc +lat_0=0 +lon_0=9 +k=1.00 +x_0=350 
+y_0=0 +ellps=bessel +units=m +nadgrids=./HeTa2010.gsb +wktext
-t_srs +proj=utm +zone=32 +ellps=GRS80 +datum=WGS84 +units=m ...

This works very well (Kudos to Frank and co-developers!!).

Now, the question came up, how to do it just the other way round, ie from
ETRS/UTM to Gauss-Krueger (DHDN) using the same grid.

Any hints?

Regards and thanks

Wolfgang

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


Re: [gdal-dev] transformation with nadgrids

2010-12-13 Thread Ralf Suhr
Hi Wolfgang,

is gdalwarp/ogr2ogr -t_srs +proj=tmerc +lat_0=0 +lon_0=9 +k=1.00 
+x_0=350 +y_0=0 +ellps=bessel +units=m +nadgrids=./HeTa2010.gsb +wktext -
s_srs +proj=utm +zone=32 +ellps=GRS80 +datum=WGS84 +units=m ... not working 
for you?

The default EPSG:25832 can also be written as +proj=tmerc +lat_0=0 +lon_0=9 
+k=0.9996 +x_0=50 +y_0=0 +ellps=WGS84 +units=m +wktext

Gr
Ralf

Am Montag 13 Dezember 2010, 12:07:51 schrieb Wolfgang Lahr:
 Hello list,
 
 I am using GDAL/OGR to transform raster and vector data from German
 Gauss-Krueger (DHDN) to ETRS/UTM like so:
 
 gdalwarp/ogr2ogr -s_srs +proj=tmerc +lat_0=0 +lon_0=9 +k=1.00
 +x_0=350 +y_0=0 +ellps=bessel +units=m +nadgrids=./HeTa2010.gsb
 +wktext -t_srs +proj=utm +zone=32 +ellps=GRS80 +datum=WGS84 +units=m
 ...
 
 This works very well (Kudos to Frank and co-developers!!).
 
 Now, the question came up, how to do it just the other way round, ie from
 ETRS/UTM to Gauss-Krueger (DHDN) using the same grid.
 
 Any hints?
 
 Regards and thanks
 
 Wolfgang
 
 ___
 gdal-dev mailing list
 gdal-dev@lists.osgeo.org
 http://lists.osgeo.org/mailman/listinfo/gdal-dev

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


[gdal-dev] gdal18dev ??

2010-12-13 Thread Agustin Lobo
Hi!

According to Synaptic (on ubuntu 10.04 64 bits) I have gdal1.7.3, but
package rgdal in R claims
GDAL runtime: GDAL 1.8dev

How can it be? I do not find any gdal18 in my system and rgdal uses
the installed gdal, as I install rgdal from source.

Thanks,

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


[gdal-dev] Question about GDALCreateGenImgProjTransformer

2010-12-13 Thread Cassanova, Bill
Hi All,

 

I am experiencing a bit of confusion regarding the use of
GDALCreateGenImgProjTransformer.

 

I am attempting to create a ProjTransformer in 2 different ways to prove
to myself that I would get the same result...The problem is I don't but
don't understand exactly why.

 

The first method I am trying is to open the GRIB2 dataset using GDALOpen
to get a GDAL dataset handle.  This handle is then passed to
GDALCreateGenImgProjTransformer as the first argument with the second
argument being passed as NULL.  The destination dataset is one created
using the well-known WGS84 Tag...For example

 

OGRSpatialReferenceH hDEST;

hDEST = OSRNewSpatialReference(NULL);

OSRSetWellKnownGeogCS(hDEST, WGS84);

OSRExportToWkt(hDEST, pszDstWKT);

OSRDestroySpatialReference(hDEST);

 

void *my_transformer = GDALCreateGenImgProjTransformer(hSrcDS,

  NULL,

  NULL,

  pszDstWKT,

  FALSE, 0.0,
0);

 

The second method passes NULL for the first argument and instead passed
the WKT instead.  The destination data set is the same as created above.

 

void *my_second_transformer = GDALCreateGenImgProjTransformer( NULL,

 
GDALGetProjectionRef( hSrcDS ),

  NULL,

 
pszDstWKT,

  FALSE,

  0.0, 0
);

 

The result of GDALGetProjectionRef( hSrcDS ) for reference :



PROJCS[unnamed,

GEOGCS[Coordinate System imported from GRIB file,

DATUM[unknown,

SPHEROID[Sphere,6371229,0]],

PRIMEM[Greenwich,0],

UNIT[degree,0.0174532925199433]],

PROJECTION[Lambert_Conformal_Conic_2SP],

PARAMETER[standard_parallel_1,25],

PARAMETER[standard_parallel_2,25],

PARAMETER[latitude_of_origin,0],

PARAMETER[central_meridian,265],

PARAMETER[false_easting,0],

PARAMETER[false_northing,0]]

 

Lambert Conformal Grid Information as extracted from gdalinfo:

Origin = (-4226106.996915472671390,7240925.337282469496131)

Pixel Size = (12191.000,-12191.000)

Corner Coordinates:

Upper Left  (-4226106.997, 7240925.337) (152d52'43.04W, 54d33'55.24N)

Lower Left  (-4226106.997, 2023177.337) (133d25'42.19W, 12d 5'14.28N)

Upper Right ( 3259167.003, 7240925.337) ( 49d15'27.05W, 57d17'54.49N)

Lower Right ( 3259167.003, 2023177.337) ( 65d 2'28.98W, 14d12'48.02N)

Center  ( -483469.997, 4632051.337) (100d30'19.44W, 40d34'21.84N)

 

 

To test both of these I pass in the known coordinates up the upper left
hand corner of the grid ( -152.877, 54.567 ) using code:

 

double x_coord = atof( argv[2] );

double y_coord = atof( argv[3] );

 

std::cout  std::fixed  original x_coord:   x_coord   y_coord:
  y_coord  std::endl;

 

result =  GDALGenImgProjTransform( my_transformer, TRUE, 1,
x_coord, y_coord,

   NULL, panSuccess);

 

  std::cout  std::fixed  calculated x_coord:   x_coord  
y_coord:   y_coord  std::endl;

 

// Reset x-coord and y-coord from command line.

  x_coord = atof( argv[2] );

  y_coord = atof( argv[3] );

 

  std::cout  std::fixed  second original x_coord:   x_coord  
y_coord:   y_coord  std::endl;

 

   esult =

  GDALGenImgProjTransform( my_second_transformer, TRUE, 1, x_coord,
y_coord,

   NULL, panSuccess);

 

   std::cout  std::fixed  second calculated x_coord:   x_coord
  y_coord:   y_coord  std::endl;

 

The results:

 

// This result is right on the money and exactly what I expect

original x_coord: -152.877000 y_coord: 54.567000

calculated x_coord: 0.016427 y_coord: -0.011915

 

// This result is not what I expect.  I would have assumed since I
created the transformer in a different way but using the WKT information
it should have been the same.

// instead I am getting the Lambert Easting and Northing in Meters.

second original x_coord: -152.877000 y_coord: 54.567000

second calculated x_coord: -4225906.734228 y_coord: 7241070.590061

 

Am I missing a step here?

 

Thanks,

Bill

 



 

William Cassanova | Senior GFS Developer | The Weather Channel |
770.226.2368 | bcassan...@weather.com

 

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

Fwd: Re: [gdal-dev] transformation with nadgrids

2010-12-13 Thread Wolfgang Lahr
Ralf,

I did not expect it to be that simple ;-). I was expecting problems
with the direction of the shifts in the HeTA-grid and therefore didn't give it 
a try.

Thanks a lot

Wolfgang



Monday, December 13, 2010, 12:47:07 PM, you wrote:


Hi Wolfgang,
is gdalwarp/ogr2ogr -t_srs +proj=tmerc +lat_0=0 +lon_0=9 +k=1.00 
+x_0=350 +y_0=0 +ellps=bessel +units=m +nadgrids=./HeTa2010.gsb +wktext 
-s_srs +proj=utm +zone=32 +ellps=GRS80 +datum=WGS84 +units=m ... not working 
for you?
The default EPSG:25832 can also be written as +proj=tmerc +lat_0=0 +lon_0=9 
+k=0.9996 +x_0=50 +y_0=0 +ellps=WGS84 +units=m +wktext
Gr
Ralf
Am Montag 13 Dezember 2010, 12:07:51 schrieb Wolfgang Lahr:
 Hello list,
 
 I am using GDAL/OGR to transform raster and vector data from German
 Gauss-Krueger (DHDN) to ETRS/UTM like so:
 
 gdalwarp/ogr2ogr -s_srs +proj=tmerc +lat_0=0 +lon_0=9 +k=1.00
 +x_0=350 +y_0=0 +ellps=bessel +units=m +nadgrids=./HeTa2010.gsb
 +wktext -t_srs +proj=utm +zone=32 +ellps=GRS80 +datum=WGS84 +units=m
 ...
 
 This works very well (Kudos to Frank and co-developers!!).
 
 Now, the question came up, how to do it just the other way round, ie from
 ETRS/UTM to Gauss-Krueger (DHDN) using the same grid.
 
 Any hints?
 
 Regards and thanks
 
 Wolfgang
 
 ___
 gdal-dev mailing list
 gdal-dev@lists.osgeo.org
 http://lists.osgeo.org/mailman/listinfo/gdal-devTitle: Re: [gdal-dev] transformation with nadgrids


Ralf,

I did not expect it to be that simple ;-). I was expecting problems with the direction of the HeTA-grid and therefore didn't give it a try.

Thanks a lot

Wolfgang



Monday, December 13, 2010, 12:47:07 PM, you wrote:





Hi Wolfgang,
is gdalwarp/ogr2ogr -t_srs "+proj=tmerc +lat_0=0 +lon_0=9 +k=1.00 +x_0=350 +y_0=0 +ellps=bessel +units=m +nadgrids=./HeTa2010.gsb +wktext" -s_srs "+proj=utm +zone=32 +ellps=GRS80 +datum=WGS84 +units=m" ... not working for you?
The default EPSG:25832 can also be written as"+proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=50 +y_0=0 +ellps=WGS84+units=m +wktext"
Gr
Ralf
Am Montag 13 Dezember 2010, 12:07:51 schrieb Wolfgang Lahr:
 Hello list,

 I am using GDAL/OGR to transform raster and vector data from German
 Gauss-Krueger (DHDN) to ETRS/UTM like so:

 gdalwarp/ogr2ogr -s_srs "+proj=tmerc +lat_0=0 +lon_0=9 +k=1.00
 +x_0=350 +y_0=0 +ellps=bessel +units=m +nadgrids=./HeTa2010.gsb
 +wktext" -t_srs "+proj=utm +zone=32 +ellps=GRS80 +datum=WGS84 +units=m"
 ...

 This works very well (Kudos to Frank and co-developers!!).

 Now, the question came up, how to do it just the other way round, ie from
 ETRS/UTM to Gauss-Krueger (DHDN) using the same grid.

 Any hints?

 Regards and thanks

 Wolfgang

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





--
mit freundlichem Gruß,
Wolfgang Lahr

Dr. Wolfgang Lahr
Büro für Geographische
Informationsverarbeitung
Gottesweg 179
50939 Köln

Tel : 0221/4600481
Fax : 0221/4759288
Mobil: 0172/6967592
Email:b...@gmx.de
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] question about ReadAsArray

2010-12-13 Thread Matt Funk
Hi,

i am a bit new to the gdal library. I am trying to process some geotiff
files in python.

1) I am making a call to the ReadAsArray data as such:
...
band = ds.GetRasterBand(1)
array = band.ReadAsArray(0,0,band.XSize,band.YSize)
print array.shape
...
which returns: (7191, 8331)
Now, what i don't understand is that when i do:
...
print 'Size is ',ds.RasterXSize,'x',ds.RasterYSize, 'x',ds.RasterCount
...
i get: Size is  8331 x 7191 x 1

So does that mean that ReadAsArray returns the an array that has its
dimensions flipped?


2) Where can i find the documentation for the ReadAsArray fcn?


thanks
matt


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


[gdal-dev] Gdalwarp -crop_to_cutline results in 0 filled raster

2010-12-13 Thread Marius Jigmond
Hi Everyone,

Has anyone experienced the problem in the subject line? I ran:

$gdalwarp -cutline shapefile -crop_to_cutline -s_srs EPSG:31700 -t_srs
EPSG:3844 source.tif dest.tif

and I ended up with a zero filled raster with the correct extent
however. Using 1.8dev. I have a rudimentary script that does the extract
by mask similar to -crop_to_cutline feature, and I was looking forward
to this much simpler approach.

-marius
-- 
Sent from Ubuntu




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


[gdal-dev] ESRI file geodatabase API finally announced

2010-12-13 Thread Bart van den Eijnden (OSGIS)
See:

http://blogs.esri.com/Dev/blogs/geodatabase/archive/2010/12/13/File-Geodatabase-API-details.aspx

Best regards,
Bart

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