1

I have GeoTIFF with these attributes:

Size is 495, 309
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6370000,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_1SP"],
    PARAMETER["latitude_of_origin",47.5],
    PARAMETER["central_meridian",4],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
PROJ.4 string is:
'+proj=lcc +lat_1=47.5 +lat_0=47.5 +lon_0=4 +k_0=1 +x_0=0 +y_0=0 +a=6370000 +b=6370000 +units=m +no_defs '
Origin = (-2963997.870569999800000,-1848004.200800000000000)
Pixel Size = (11975.754948303031000,11961.191901488673000)
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-2963997.871,-1848004.201) ( 24d36'23.04"W, 26d22' 5.88"N)
Lower Left  (-2963997.871, 1848004.097) ( 45d39'38.75"W, 55d16'35.82"N)
Upper Right ( 2964000.829,-1848004.201) ( 32d36'23.13"E, 26d22' 5.85"N)
Lower Right ( 2964000.829, 1848004.097) ( 53d39'38.88"E, 55d16'35.76"N)
Center      (   1.4791350,  -0.0520200) (  4d 0' 0.07"E, 47d30' 0.00"N)
Band 1 Block=495x4 Type=Float32, ColorInterp=Gray
  NoData Value=9.969209968386869e+036

and I want to overlay Natural Earth shapefile.

So I have to crop and reproject global Natural Earth shapefile to this GeoTIFF sample, however I can't understand how to do this with such projection and was wondering if I can instruct somehow ogr2ogr to use my GeoTIFF as reference to crop and reproject the shapefile?

Also I'm interested in other approaches if ogr2ogr does not have this feature.

2

If the Proj4 string that is stored into the GeoTIFF tags is correct this should be an easy case. Without cropping it would look like

ogr2ogr -f "ESRI Shapefile" -s_srs epsg:4326 -t_srs "+proj=lcc +lat_1=47.5 +lat_0=47.5 +lon_0=4 +k_0=1 +x_0=0 +y_0=0 +a=6370000 +b=6370000 +units=m +no_defs" projected_out.shp natural_earth_in.shp

For cropping to the destination extents, you can see the corner coordinates in the gdalinfo report. Use the in ogr2ogr parameter -clipdst xmin ymin xmax ymax. Read more from the manual http://www.gdal.org/ogr2ogr.html.

EDIT

Full command which reprojects and clips countries of the world to target extents and target projection is

ogr2ogr -f "ESRI Shapefile" -s_srs epsg:4326 -t_srs "+proj=lcc +lat_1=47.5 +lat_0=47.5 +lon_0=4 +k_0=1 +x_0=0
 +y_0=0 +a=6370000 +b=6370000 +units=m +no_defs" projected_out.shp ne_110m_admin_0_countries.shp -clipdst -2963997 18480
04 2964000 -1848004 --config OGR_ENABLE_PARTIAL_REPROJECTION TRUE -skipfailures

Notice the use of GDAL config option OGR_ENABLE_PARTIAL_REPROJECTION which is compulsory in this case because it is impossible to reproject the whole source dataset into target projection which is valid only for a limited area. Also option -skipfailures was needed either because source data has topology errors or because self-intersections were created in the reprojection process. Result was still not perfect - Russia is missing and Antarctica is covering the whole Europe which means that some further investigations would be needed. That is worth another question, though. The light background in the image below is the miss-projected Antarctica.

Directly clipped result

  • Projecting whole globe to destination projection is not possible, so I have to crop the source first, and that is where I stopped in the first place as corner coordinates are not forming rectangle in Natural Earth shapefile source projection. Any suggestions? – zetah Aug 27 '14 at 17:37
  • I can have a try. Which Natural Earth shapefile? – user30184 Aug 27 '14 at 18:09
  • Thanks, it was countries shapefile, but I managed to get correct output by first cropping to outer rectangle of GeoTIFF file, then reprojecting and then clipping to projected coordinates. I was just hoping that straightforward approach of using georeferenced file as reference exists in gdal/org system. – zetah Aug 27 '14 at 19:33
  • Clipping to extents of a raster file is not supported directly but clipping to vector layer is. You should just convert the extents of your raster layer into a polygon by digitizing a bounding polygon or by creating BBOX for example with gdaltindex utility. GDAL is quite capable. I edited my answer to show how to avoid the first cropping that you did. – user30184 Aug 27 '14 at 20:34
  • Nice tip. I tried to reproject and clip at the same time, but I put -clipdst before the filenames, so it didn't work, and I didn't read the error message carefully about partial reprojection possiblity with OGR_ENABLE_PARTIAL_REPROJECTION nor I would have known how to define that option ;) – zetah Aug 27 '14 at 22:10

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.