3

I would like to know correct syntax of R {gdalUtils} package (version 2.0.1.7) to georeference Shapefiles using corner points as GCPs. (These Shapefiles have y=height and x=width information, without datum / projection.)

On OSGeo4W shell, ogr2ogr works fine and it returns georeferenced Shapefile.

ogr2ogr -gcp 0 0 40 40 -gcp 0.263693 0 145 40 -gcp 0.263693 -0.237619 145 -50 -gcp 0 -0.237619 40 -50 -f "ESRI shapefile" Myreferenced.shp My.shp

Trying to do the same by gdalUtils::ogr2ogr(),

library(gralUtils)
src_datasource_name # input Shapefile
dst_datasource_name # output Shapefile (georeferenced)

# ogrinfo() works fine
ogrinfo(datasource_name = src_datasource_name, al = TRUE, so = TRUE)

# ogr2ogr() fails
ogr2ogr(src_datasource_name= src_datasource_name,
        dst_datasource_name= dst_datasource_name,
        f= "ESRI Shapefile",
        gcp = list(1:4,
                   c(0, 0, 40, 40),
                   c(0.263693, 0, 145, 40),
                   c(0.263693, -0.237619, 145, -50),
                   c( 0, -0.237619, 40, -50)
                  )
       )

What I am struggling in the ogr2ogr() syntax is how to provide multiple GCP points. As you see above, I tried to give GCPs as a list, which did not work.

Does anybody know how to set GCPs with 4 corner points in gdalUtils? Unfortunately use cases of gdalUtils are relatively scarce on internet.


About my Shapefile Data

  • They were originally hand-drawn by MS-Powerpoint and saved as .eps file. It has bounding box with known LatLong corner points (UL, UR, LR, and LL). The .eps file was then converted to .shp file.
  • Conversion from .eps to .shp changes everything into Line feature (Multilines Shapefile).
  • This Shapefile does not have projection information but its dimension can be read by ogrinfo. It is (0, 0) at the UL corner, with height (eg y= -0. 238) and width (x= 0.263). Measured unit is meter. These y and x values can vary among drawings.

Error message from above ogr2ogr() is

[1] "ERROR 1: Failed to compute GCP transform: Not enough points available"

However increasing points does not seem a solution, as ogr2ogr on OSGeo4W shell works fine on the same data and GCP sets. This error is due to my wrong code, but I do not know how to fix it.


[EDIT] by QGIS Raster georeferencer, .points would look like this image.

gcpimage

mapX,mapY,pixelX,pixelY,enable
40,40,0,0,1
145,40,0.26300000000000001,0,1
40,-50,0,-0.23799999999999999,1
145,-50,0.26300000000000001,-0.23799999999999999,1

[EDIT2] Another way of applying GCPs:

ogr2ogr(src_datasource_name= src_datasource_name,
        dst_datasource_name= dst_datasource_name,
        f= "ESRI Shapefile",
        gcp = c(0,                0,  40,  40),
        gcp = c(0.263693,         0, 140,  40),
        gcp = c(0.263693, -0.237619, 140, -50),
        gcp = c(0,        -0.237619,  40, -50)
        )

This way R returned following error

Error in ogr2ogr(src_datasource_name = src_datasource_name, dst_datasource_name = dst_datasource_name, : formal argument "gcp" matched by multiple actual arguments

I don't know what it means...


[EDIT3] Suggested approach to put variables into -optfile.

ogr2ogr(src_datasource_name= src_datasource_name,
        dst_datasource_name= dst_datasource_name,
        f= "ESRI Shapefile",
        optfile = "C:/temp/my_gcp.txt"
        )

Below is my tested formats (my_gcp.txt) yet successful...

# trial 1
"-gcp 0 0 40 40 -gcp 0.263693 0 145 40 -gcp 0.263693 -0.237619 145 -50 -gcp 0 -0.237619 40 -50"
# trial 2
"gcp = c(0, 0, 40, 40,), gcp = c(0.263693, 0, 145, 40), gcp = c(0.263693, -0.237619, 145, -50), gcp = c(0, -0.237619, 40, -50)"

Error message was:

Error in ogr2ogr(src_datasource_name = src_datasource_name, dst_datasource_name = dst_datasource_name,  :
unused argument (optfile = "C:/temp/my_gcp.txt")
  • Print out the "gcp" once you have set the value with your list. How does it look like? – user30184 Nov 4 '16 at 11:34
  • @user30184, Thank you for your response. I was not so sure what would be the best to print out to see its look, so I used QGIS to try to make somewhat similar. Maybe I am mistaken? – Kazuhito Nov 4 '16 at 12:16
  • Make your code to print the contents of the variable. It may be -gcp followed by many mapX,mapY,pixelX,pixelY, while you will need many times "-gcp mapX,mapY,pixelX,pixelY". – user30184 Nov 4 '16 at 13:06
  • Got that, but the suggested way, putting -gcp (or gcp without hyphen) commands in row (which is, genuine GDAL style) was what I had tried first. It seemed enumeration of objects was not acceptable to R. – Kazuhito Nov 4 '16 at 13:19
  • @user30184 Apologies. Forgot to add the error message I got with multiple gcp lines. Just added as [EDIT 2] above. Hope this clarifies the situation better... – Kazuhito Nov 4 '16 at 14:18

Your Answer

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

Browse other questions tagged or ask your own question.