[GIS] GDAL GetGeoTransform returning lat/lng instead of affine co-ordinates

ccoordinate systemgdal

I'm using the GDAL libraries via the c# wrapper and have hit a problem when using geotiff's where the co-ordinates are from a state plane zone. Previously I was having trouble getting the correct affine co-ordinates converted to lat/lng.

Having figured out the problem there, I now find that when I get the GeoTransform values from the dataset, it contains lat/lng rather than the affine co-ordinates. This makes it harder to calculate he co-ordinates of each corner of a raster. (In fact, I'm not sure how to do that).

On top of that, I have to somehow detect that this raster is going to give me lat/lng rather than affine co-ordinates. The API doc doesn't mention that lat/lng may turn up in the GetGeoTransform()…should it?

My code:

double[] geot = new double[6];
ds.GetGeoTransform(geot);

returns:
-77.5286920001666,
1,
0,
38.728995261289775,
0,
-1

Is there some way I can either find the affine co-ordinates using the lat/lng, or even better, get the dataset object to give me the affine co-ordinates instead of the lat/lng?

EDIT: Adding GetProjectionRef() here:

PROJCS["NAD_1983_HARN_StatePlane_Virginia_North_FIPS_4501_Feet",
GEOGCS["NAD83(HARN)",
DATUM["NAD83_High_Accuracy_Reference_Network",
SPHEROID["GRS 1980",6378137,298.2572221010002,
AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","6152"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4152"]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",39.2],
PARAMETER["standard_parallel_2",38.03333333333333],
PARAMETER["latitude_of_origin",37.66666666666666],
PARAMETER["central_meridian",-78.5],
PARAMETER["false_easting",37673535.76388889],
PARAMETER["false_northing",21527734.72222222],
UNIT["us_survey_feet",0.3048006096012192],
AUTHORITY["EPSG","2853"]]

Best Answer

Andrew,

I think you're thinking wrong and it may get you into trouble. You need to stop thinking in terms of your dataset and start thinking in terms of GDAL and GIS. Your dataset data is in some projection, in this case State Plane. That means you use state plane (ft here) coordinates to access the data. The transform (GetGeoTransform) converts between your dataset projection coordinates and the pixel indices into the raster 2D array of values. You can call OSGeo.GDAL.Gdal.InvGeoTransform() to invert this transform to convert pixels to your dataset projection.

Now, if you want lat/lon, you need to pick a geodetic coordinate system (like WGS84). Then you construct the SRS for that system (using the WKT string), as well as get the SRS for your dataset (GetProjectionRef). Once you have those, you are in fat city.

You construct the transform to go from geodetic (wgs84) -> projected (state plane) using OGRCreateCoordinateTransformation, and you use that to find the state plane values for your lat/lon. You can construct the inverse transform if you need to go from lat/lon to state plane.

Once you have state plane, you use the affine to go from state plane to pixels. And you can use the inverse affine to go from pixels to state plane. So, if your raster is (say) 4000 rows by 1200 columns, you can use the corners (0, 1199 and 3999) into the inverse transform to get the state plane, and then use the transform to get them back to lat/lon.

The point is, you stop worrying about your dataset SRS and the values in the affine projection. Always think in terms of the SRS that you want to work in (e.g. WGS84), your dataset's SRS (which you get from GetProjectionRef), and then pixels (which you get from the affine transform, GetGeoTransform). Just convert between those and use them where you need them.

Clear as mud?

-reilly.

Related Question