1

I am very new to gis, and while most of the stuff I have is working well, this thing has me stumped.

I'm using GDAL libraries and the c# wrapper to access various gis functions. I'm processing ortho images to merge a few tiles, convert to WGS84 and then I chop the merged image into a bunch of small tiles. For the most part, this is working well, but I've come across this image which I can't get the right co-ordinate for.

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"]]

I use the following code to extract co-ordinates (WGS84_STR is the WGS84 WKT):

string proj = ds.GetProjectionRef();
SpatialReference old_cs = new SpatialReference(proj);
SpatialReference new_cs = new SpatialReference(WGS84_STR);
trans = new CoordinateTransformation(old_cs, new_cs);

...then a method to return the coordinates:

...
ds.GetGeoTransform(adfGeoTransform);
dfGeoX = adfGeoTransform[0] + adfGeoTransform[1] * x + adfGeoTransform[2] * y;
dfGeoY = adfGeoTransform[3] + adfGeoTransform[4] * x + adfGeoTransform[5] * y;
trans.TransformPoint(adfGeoTransform,  dfGeoX, dfGeoY, 0);
...

where x and y are the image corner offsets.

I get lat of: -14.921860233519684, and lng of: -130.08927419763472 ...which is in the Pacific Ocean somewhere, not Northern Virginia. The tile is also rotated. I wonder if someone could point out where I've gone wrong?


Source of problem found It turns out the GEOTIFF_CSV environmental variable was not set. In c#, the following code: OSGeo.GDAL.Gdal.SetConfigOption("GEOTIFF_CSV", Path.GetDirectoryName(Application.ExecutablePath) + @"\gdal-data"); makes a very big difference. Now the geoTransform contains the actual lat/long.

  • What is the extent of the input raster? – mkennedy Mar 16 '15 at 13:50
  • The tile is 5000x5000 pixels, 1 foot resolution. the actual bounding coordinates are WEST LONGITUDE = - 78.66761249 º W NORTH LATITUDE =39.33129819 º N EAST LONGITUDE = - 76.51363560 º W SOUTH LATITUDE =37.76418697 º N – Andrew C Mar 16 '15 at 14:34
  • May be related to this GDAL ticket trac.osgeo.org/gdal/ticket/4954 – user30184 Mar 16 '15 at 14:58
  • If that's the actual extent (and not just the lat/lon extent), it means the raster isn't in a state plane zone, but is in lat/lon. – mkennedy Mar 16 '15 at 16:09
  • 1
    Just discovered that fwtools - tiffinfo reports a state plane zone of 4501, and when I plug this info into link , I get the right co-ordinates. Soooo, I used the API SetStatePlane to set the plane info, and the measurement info, and bingo - right result. It seems the gdal library can't detect the plane zone by itself in this case. (there is no GetStatePlane either). – Andrew C Mar 18 '15 at 16:31
1

Source of problem found. It turns out the GEOTIFF_CSV environmental variable was not set. In c#, the following code:

OSGeo.GDAL.Gdal.SetConfigOption("GEOTIFF_CSV", Path.GetDirectoryName(Application.ExecutablePath) + @"\gdal-data"); 

makes a very big difference. Now the geoTransform contains the actual lat/long.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.