4

I'm having real problems trying to fit a GeoPDF to a set of shapefiles using gdaltransform and gdalwarp. The problem I can't seem to fix is that each time I run the commands, whether in QGIS or command line, the new GeoTiff is always off by a few hundred meters. I've read about using the +towgs84 options when setting the CRS, but the online documentation for this option set is very vague about how to set the 3-7 parameters and which place to put them.

In complete and utter despondency, I started to think maybe I'm not using the correct CRS in the first place. When running gdalinfo, I see the multiple EPSG codes. I had been using EPSG:4267, thinking this was the actual CRS of the GeoPDF. But if EPSG:4267 is indeed the correct CRS, then why does the file keep showing an offset of ~188 meters when trying to match it to my shapefiles in EPSG:4326?

Here's the gdalinfo output of my GeoPDF file. Can someone please identify what is the correct CRS of the file, and how I can use the CRS to eliminate the offset?

Driver: PDF/Geospatial PDF Files: ~/Desktop/myfile.pdf Size is 4470, 3375 Coordinate System is: PROJCS["UTM Zone 18, Northern Hemisphere", GEOGCS["NAD27", DATUM["North_American_Datum_1927", SPHEROID["Clarke 1866",6378206.4,294.9786982138982, AUTHORITY["EPSG","7008"]], AUTHORITY["EPSG","6267"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4267"]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",-75], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["Meter",1]] GeoTransform = 786521.192808, 8.469778080000001, 0.15067296 2179199.881801, 0.14403264, -8.469838560000001 GCP Projection = PROJCS["UTM Zone 18, Northern Hemisphere", GEOGCS["NAD27", DATUM["North_American_Datum_1927", SPHEROID["Clarke 1866",6378206.4,294.9786982138982, AUTHORITY["EPSG","7008"]], AUTHORITY["EPSG","6267"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4267"]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",-75], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["Meter",1]] GCP[ 0]: Id=1, Info= (207.8125,2466.92708333333) -> (788652.79799,2158337.734581,0) GCP[ 1]: Id=2, Info= (211.197916666667,287.760416666667) -> (788355.62363,2176795.854646,0) GCP[ 2]: Id=3, Info= (3308.59375,288.020833333333) -> (814587.658583,2177239.223702,0) GCP[ 3]: Id=4, Info= (3307.8125,2467.70833333333) -> (814911.952376,2158777.947799,0) GCP[ 4]: Id=5, Info= (2275.78125,1377.34375) -> (806001.09572,2167857.005616,0) GCP[ 5]: Id=6, Info= (1242.70833333333,1377.08333333333) -> (797252.559747,2167709.741736,0) Metadata: AUTHOR=NGA CREATION_DATE=D:20071008152059-04'00' KEYWORDS=003 NEATLINE=POLYGON ((788648.68490375 2158330.85305763,788592.08154525 2161760.771554,788380.2253165 2174662.13561613,788346.57604175 2176801.72714438,802820.27008 2177043.44654713,812895.73768325 2177216.99106225,814591.898970625 2177245.83509875,814733.508010375 2169285.51026463,814856.815543875 2162105.92365575,814916.143021875 2158770.92472275,802630.671176375 2158555.38531638,793697.623632125 2158407.88726563,789672.19489725 2158343.84562738,788648.68490375 2158330.85305763,788648.68490375 2158330.85305763,788648.68490375 2158330.85305763,788648.68490375 2158330.85305763)) TITLE=Map Edition Corner Coordinates: Upper Left ( 786521.193, 2179199.882) ( 72d16' 1.60"W, 19d41'19.09"N) Lower Left ( 787029.714, 2150614.177) ( 72d15'59.83"W, 19d25'49.84"N) Upper Right ( 824381.101, 2179843.708) ( 71d54'22.52"W, 19d41'18.91"N) Lower Right ( 824889.622, 2151258.003) ( 71d54'22.82"W, 19d25'49.93"N) Center ( 805705.407, 2165228.942) ( 72d 5'11.67"W, 19d33'34.77"N) Band 1 Block=4470x1 Type=Byte, ColorInterp=Red Band 2 Block=4470x1 Type=Byte, ColorInterp=Green Band 3 Block=4470x1 Type=Byte, ColorInterp=Blue

0
7

Your data seems to be in Haiti, using the North American Datum of 1927.

QGIS and GDAL use the datum grid shift files for NAD27 for USA and Canada, but unfortunately your data is just outside the conus grid file. hence all your reprojections get a null shift, which is definitely wrong.

Clifford J. Mugnier published a Grids & Datums article on Haiti in September 2011, giving these shift values:

from NAD1927 to WGS84: DeltaX =  –3 m ± 3 m, DeltaY = +142 m ± 9 m, 
DeltaZ = +183 m ± 12 m. Curiously, the shift parameters that work out 
for  the single Haïtian point, “FORT NATIONAL”  from NAD1927  to 
WGS84 are: DeltaX =  –13 m, DeltaY =  –95 m, DeltaZ =  –197 m. 

Your data is in NAD27/UTM 18N, but EPSG:26718 fails for the reasons given above, and you have to build -s_srs and -t_srs for every GDAL-based operation manually:

+proj=utm +zone=18 +ellps=clrk66 +towgs84=-3,142,183,0,0,0 +units=m +no_defs
+proj=utm +zone=18 +ellps=clrk66 +towgs84=-13,-95,-197,0,0,0 +units=m +no_defs

Taking the Trou du Nord map from http://www.lib.utexas.edu/maps/topo/haiti/ and the first datum shift, it fits quite good to an OpenStreetMap background, with the neatline in red:

enter image description here

2
  • BRAVO! Let's say I want to tweak the +towgs84 values for this particular map -- let's say 24 additional meters to the west and 50 meters to the south. How would I do this using the 1st data shift you provided? Sep 19 '17 at 18:28
  • The shift values count from the Earth`s center, so no direct calculation is possible. Just use try&error. add 20 meters to DeltaX, and see how the image moves in the right or wrong direction. Change the value until you are close to where you want, then do the same with deltaY, then with deltaZ. If the map gets offset in other parts of the image, you might add the rotation parameters in the same way.
    – AndreJ
    Sep 20 '17 at 5:41
0

Pasting the WKT into http://prj2epsg.org gives the plausible sounding EPSG:2032

PROJCS["NAD27(CGQ77) / UTM zone 18N", 
  GEOGCS["NAD27(CGQ77)", 
    DATUM["North American Datum 1927 (CGQ77)", 
      SPHEROID["Clarke 1866", 6378206.4, 294.9786982138982, AUTHORITY["EPSG","7008"]], 
      AUTHORITY["EPSG","6609"]], 
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], 
    UNIT["degree", 0.017453292519943295], 
    AXIS["Geodetic longitude", EAST], 
    AXIS["Geodetic latitude", NORTH], 
    AUTHORITY["EPSG","4609"]], 
  PROJECTION["Transverse Mercator", AUTHORITY["EPSG","9807"]], 
  PARAMETER["central_meridian", -75.0], 
  PARAMETER["latitude_of_origin", 0.0], 
  PARAMETER["scale_factor", 0.9996], 
  PARAMETER["false_easting", 500000.0], 
  PARAMETER["false_northing", 0.0], 
  UNIT["m", 1.0], 
  AXIS["Easting", EAST], 
  AXIS["Northing", NORTH], 
  AUTHORITY["EPSG","2032"]]
0

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.