3

I have one elevation file (SRTM) and one satellite imagery file (sentinel-2), both are in GeoTIFF format. The elevation file has a much larger area covered. Here is the gdalinfo of the elavation file:

> Driver: GTiff/GeoTIFF Files: C:\DEV\3D\Assets\Heightmaps\nv\full.tif
> Size is 10801, 7201 Coordinate System is: GEOGCS["WGS 84",
>     DATUM["WGS_1984",
>         SPHEROID["WGS 84",6378137,298.257223563,
>             AUTHORITY["EPSG","7030"]],
>         AUTHORITY["EPSG","6326"]],
>     PRIMEM["Greenwich",0],
>     UNIT["degree",0.0174532925199433],
>     AUTHORITY["EPSG","4326"]] Origin = (-2.000138888888889,36.000138888888884) 
> Pixel Size =(0.000277777777778,-0.000277777777778) 
> Metadata:   AREA_OR_POINT=Area
> Image Structure Metadata:   INTERLEAVE=BAND Corner Coordinates: 
> Upper Left  (  -2.0001389,  36.0001389) (  2d 0' 0.50"W, 36d 0' 0.50"N)
> Lower Left  (  -2.0001389,  33.9998611) (  2d 0' 0.50"W,33d59'59.50"N) 
> Upper Right (   1.0001389,  36.0001389) (  1d 0'0.50"E, 36d 0' 0.50"N) 
> Lower Right (   1.0001389,  33.9998611) (  1d 0' 0.50"E,33d59'59.50"N) 
> Center(  -0.5000000,  35.0000000) ( 0d30' 0.00"W, 35d 0' 0.00"N) 
> Band 1 Block=10801x1 Type=Int16,ColorInterp=Gray

and the satellite image:

>   Driver: GTiff/GeoTIFF Files:
> C:\DEV\3D\Assets\Imagery\new\sub_big_merge.tif Size is 15017, 8131
> Coordinate System is: PROJCS["WGS 84 / UTM zone 30N",
>     GEOGCS["WGS 84",
>         DATUM["WGS_1984",
>             SPHEROID["WGS 84",6378137,298.257223563,
>                 AUTHORITY["EPSG","7030"]],
>             AUTHORITY["EPSG","6326"]],
>         PRIMEM["Greenwich",0,
>             AUTHORITY["EPSG","8901"]],
>         UNIT["degree",0.0174532925199433,
>             AUTHORITY["EPSG","9122"]],
>         AUTHORITY["EPSG","4326"]],
>     PROJECTION["Transverse_Mercator"],
>     PARAMETER["latitude_of_origin",0],
>     PARAMETER["central_meridian",-3],
>     PARAMETER["scale_factor",0.9996],
>     PARAMETER["false_easting",500000],
>     PARAMETER["false_northing",0],
>     UNIT["metre",1,
>         AUTHORITY["EPSG","9001"]],
>     AXIS["Easting",EAST],
>     AXIS["Northing",NORTH],
>     AUTHORITY["EPSG","32630"]] 
> Origin = (614090.000000000000000,3987110.000000000000000) 
> Pixel Size =(10.000000000000000,-10.000000000000000) 
> Metadata:  
> AREA_OR_POINT=Area Image Structure Metadata:   INTERLEAVE=PIXEL Corner
> Coordinates: 
> Upper Left  (  614090.000, 3987110.000) (1d44' 1.72"W,36d 1'18.55"N) 
> Lower Left  (  614090.000, 3905800.000) (1d44'43.19"W, 35d17'19.93"N) 
> Upper Right (  764260.000, 3987110.000) (0d 4' 6.98"W, 35d59'33.60"N) 
> Lower Right (  764260.000, 3905800.000) (0d 5'42.84"W, 35d15'37.76"N) 
> Center      (  689175.000, 3946455.000) (0d54'37.94"W, 35d38'37.76"N) 
> Band 1 Block=15017x1 Type=Byte,
> ColorInterp=Red Band 2 Block=15017x1 Type=Byte, ColorInterp=Green Band
> 3 Block=15017x1 Type=Byte, ColorInterp=Blue

Now, how do I get a subset of the elevation file that exactly matches the same area of the satellite image with gdal_translate or any other tool (since as you can see the two files have differents projections, resolutions, units ... etc)?

In the meantime I have tried these two "solutions"

Use GDAL command line to copy projections

Use gdalsrsinfo to get the srs of the tiff that still has the projection:

gdalsrsinfo -o wkt tiffwithsrs.tiff

Then copy the output and use gdal_translate to apply it to a new tiff:

gdal_translate -a_srs '...' tiffwithoutsrs.tif newfixedtif.tif

just paste your projection after the -a_srs option wich did nothing .

and this which completely broke my file.

5

First of all, I'd move the two raster layers in the same working directory in order to skip the long paths in command line.

We can transform the elevation layer to the same CRS of the satellite image, i.e. EPSG:32630:

gdalwarp -t_srs EPSG:32630 full.tif full_32630.tif

Then we can compute the common area between the two rasters:

gdaltindex elevation_tile.shp full_32630.tif
gdaltindex sat_image_tile.shp sub_big_merge.tif

We can write an OGR VRT layer called e.g. tiles.vrt like this:

<OGRVRTDataSource>
    <OGRVRTLayer name="elevation_tile">
        <SrcDataSource>elevation_tile.shp</SrcDataSource>
    </OGRVRTLayer>
    <OGRVRTLayer name="sat_image_tile">
        <SrcDataSource>sat_image_tile.shp</SrcDataSource>
    </OGRVRTLayer>
</OGRVRTDataSource>

and calculate the intersection:

ogr2ogr intersection.shp tiles.vrt -dialect sqlite -sql "SELECT ST_Intersection(elevation_tile.geometry, sat_image_tile.geometry) AS geometry FROM elevation_tile, sat_image_tile"

Finally, we can crop the two raster layers with intersection.shp and set the target resolution as 30 m (because 0.000277777777778° =~ 30m is the lowest resolution, it should be the target one):

gdalwarp -cutline intersection.shp -tr 30 30 full_32630.tif elevation_cropped_30m.tif
gdalwarp -cutline intersection.shp -tr 30 30 sub_big_merge.tif sat_image_cropped_30m.tif

Result: elevation_cropped_30m.tif and sat_image_cropped_30m.tif spatial extent and resolution should match as requested.

  • woow, id didn't expected that complexity XD, i have a naive idea : the satellite image bounds are (614090.000, 3987110.000 764260.000, 3905800.000) isn't possible to simply convert that coordinates to the projection of the elevation file, then use gdal_translate on the elevation file with -projwin to get a subset that much the image ? – user2475096 Apr 3 '17 at 19:44
  • @user2475096 There are different ways to achieve the result. You can crop your data manually, but misalignments are behind the corner... Try it and then let me know if it works fine! – Antonio Falciano Apr 4 '17 at 8:20
1

If you can use Python and its GDAL bindings, you can have a look at my reproject_image_to_raster method. Basically, it takes one image as the master (in this case, your S2), and it would reproject the slave to match extension and output resolution. I think you only need access to the gdal (e.g. you need an import gdal).

  • i don't want to match the resolution, only bounds or area if you want – user2475096 Apr 7 '17 at 1:01
  • The slave final resolution is an option. By default, it makes it match the input (I think...), so just set that option to the original raster resolution. – Jose Apr 7 '17 at 17:05

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.