2

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"]]
  • What is the output of ds.GetProjectionRef()? – user2856 Mar 21 '15 at 10:50
  • Hi Luke, it was in the"previously" link, but Ill stick in here too. – Andrew C Mar 21 '15 at 15:04
  • One thing that strikes me: you mention "affine coordinates" several times. An "affine* is a type of transformation, a transformation between coordinate systems, not a type of coordinate. Not claiming that is going to help you with an answer -- just an observation. – Martin F Mar 23 '15 at 0:46
  • And a belated welcome to GIS SE :-) – Martin F Mar 23 '15 at 0:57
  • Thanks @Martin F. You are right, and it's because I don't really know the right terminology :P The API for GetGeoTransform() says "Fetch the affine transformation coefficients", so I just went with it. I was expecting XY coordinates, and got lat/lng. And thanks for the welcome :) – Andrew C Mar 24 '15 at 13:02
2

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.

  • Hi reilly. I provided a "previous" link which may help explain the issue. The API documentation says GetGeoTransform returns affine co-ordinates, however this is not true for this raster for some reason. So when you say "The transform (GetGeoTransform) converts between your dataset projection coordinates and the pixel indices into the raster 2D array of values", I wish this was true, but as you can see in my post, it doesn't. 38.728995261289775, -77.5286920001666 is the lat/lng of the TL corner of the raster. GetGeoTransform normally returns affine co-ordinates, which are easy to use. – Andrew C Mar 22 '15 at 17:08
  • And what I'm saying is, that is how it works. If it isn't working this way in your case, then either your GeoTiff projection is wrong, or else your GDAL environment is wrong. Have you used the command line tools to interrogate your GeoTiff? You can use gdalinfo and gdaltransform to test your dataset to see if it is incorrectly formed. Plus, running the tools from the command line will confirm for you whether your environment is set up properly. – Reilly Mar 23 '15 at 16:56
  • Good point @Reilly. Using the following command in fwtools (gdalwarp do_n16_7964_20.2013.tif 10seh475307_wgs84.tif -t_srs "+proj=longlat +ellps=WGS84"), the raster gets rotated and placed in the Pacific Ocean. Same result I was getting initially (even though the fwtools env is separate to the dev env). But I do get the feeling I'm missing something... – Andrew C Mar 24 '15 at 2:10
  • I have never used the fwtools and never saw a reason. I just get the latest gdal release. It contains all the tools and the C# wrappers built for AnyCPU. I was kind of uncertain as to how often fwtools was updated. In any case, that means something is amiss with your tiff. Another trick is to load it into QuantumGIS. That is free and uses GDAL under-the-hood. I'm fairly certain that you'll see the problem persisting. At that point, you have to decide how to handle this bad dataset. I fear that if you code around it, you'll be sorry later! – Reilly Mar 25 '15 at 12:03
  • Another thing, the GDAL distro has a batch file that opens a command prompt with all the paths set up correctly. That batch file saved my bacon, because things weren't working in my code, but were on the command line. By looking at the batch file, I realized that I wasn't setting up all the required GDAL environment variables. Can you tell I'm a total GDAL fanboy? – Reilly Mar 25 '15 at 12:08

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.