trying to get lat/lon out of sdts dem

[nakiya123 at y...]

--- In gdal-dev at y..., Frank Warmerdam <warmerdam at p...> wrote:
> nakiya123 at y... wrote:
> 
> > Howdy,
> > 
> > I am trying to extract the lat/lon coordinates out of sdts dem's. I
> > have looked into and tried to used the OGRSpatialReference class but
> > have not seen any clear way to do this. How is this done?
> > 
> > Thanks in advance
> 
> 
> Hi,
> 
> Are the DEMs in lat/long already, or UTM? If they are in lat/long you
> can just ignore the OGRSpatialReference. Determine the georeferenced
> coordinate of the top left corner of a pixel as follows:
> 
> 
> dfGeorefX = padfGeoTransform[0] + padfGeoTransform[1] * nPixel
> + padfGeoTransform[2] * nLine;
> dfGeorefY = padfGeoTransform[3] + padfGeoTransform[4] * nPixel;
> + padfGeoTransform[5] * nLine;
> 
> If the DEM is in UTM, the above will give you the corner in UTM
coordinates.
> You can then translate this to LatLong using something like the
following
> to create a reprojector:
> 
> #include "gdal_priv.h"
> #include "ogr_spatialref.h"
> 
> OGRSpatialReference *poUTM = new OGRSpatialReference();
> const char *pszUTM_WKT = poDS->GetProjectionRef();
> 
> poUTM->importFromWkt( (char **) &pszUTM_WKT );
> 
> OGRSpatialReference *poLL = poUTM->CloneGeogCS();
> 
> OGRCoordinateTransform *poTransform;
> 
> poTransform = OGRCreateCoordinateTransformation( poUTM, poLL );
> 
> 
> Then invoke it on points like this:
> 
> 
> poTransform->Transform( 1, &dfGeorefX, &dfGeoref );
> 
> I hope this helps. Let me know if you need more details.
> 
> Best regards,
> 
> --

!THANKS FOR THE REPLY!

Just a quick note, poTransform->Transform should be
poTransform->Transformation I think. That worked for me.

I had already tried this and decided I was doing something wrong
because the values of the dfGeoref's. 
i.e. I am using a USGS SDTS DEM in UTM format and after displaying the
results of this code I get:

dfGeorefX = -115.489
dfGeorefY = 0

my GetProjectRef is:
PROJCS["UTM Zone 12, Northern
Hemisphere",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Clarke
1866",6378206.4,294.978698213898,AUTHORITY["EPSG",7008]],TOWGS84[-3,142,183,0,0,0,0],AUTHORITY["EPSG",6267]],PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG",9108]],AXIS["Lat","NORTH"],AXIS["Long","EAST"],AUTHORITY["EPSG",4267]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0]]

and the actual lat/lon of the image is in Arizona with a lat/lon of
aprox 33.xxx/111.xxx.

I'm not clear as to what the results of the Transform function are?
I thought they may be the bottom left corner of the UTM Zone but
that's not very useful unless I'm missing something?

Here's the code I used.


double dfGeorefX=0.0, dfGeorefY=0.0;

OGRSpatialReference *poUTM = new OGRSpatialReference();
const char *pszUTM_WKT = poDataset->GetProjectionRef();

cout << endl << pszUTM_WKT << endl << endl;

poUTM->importFromWkt((char**) &pszUTM_WKT);
OGRSpatialReference *poLL = poUTM->CloneGeogCS();

OGRCoordinateTransformation *poTransform;

poTransform = OGRCreateCoordinateTransformation(poUTM, poLL);

poTransform->Transform(1, &dfGeorefX, &dfGeorefY);//dfGeorefY instead

cout << endl << dfGeorefX << " " << dfGeorefY << endl << endl;

thanks again!!!
---------------------------------------+--------------------------------------
> I set the clouds in motion - turn up | Frank Warmerdam, warmerdam at p...
> light and sound - activate the windows | http://pobox.com/~warmerdam
> and watch the world go round - Rush | Geospatial Programmer for Rent