3

I have downloaded tiles from http://srtm.csi.cgiar.org/srtmdata/

http://srtm.csi.cgiar.org/srtmdata/ and unpacked the *.zip files into a directory. Then I ran, as has been suggested here: CGIAR SRTM to xyz using GRASS GIS, with sampling, the commands

#create mosaik (optionally reproject on the fly with -t_srs)
gdalwarp *.tif mosaik.tif
gdalbuildvrt srtm.vft mosaic.tif

After this, I can see the following:

Driver: VRT/Virtual Raster
Files: srtm.vft
       mosaic.tif
Size is 48001, 24001
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 = (-10.000416206603404,55.000417150924136)
Pixel Size = (0.000833333333333,-0.000833333333333)
Corner Coordinates:
Upper Left  ( -10.0004162,  55.0004172) ( 10d 0' 1.50"W, 55d 0' 1.50"N)
Lower Left  ( -10.0004162,  34.9995838) ( 10d 0' 1.50"W, 34d59'58.50"N)
Upper Right (  30.0004171,  55.0004172) ( 30d 0' 1.50"E, 55d 0' 1.50"N)
Lower Right (  30.0004171,  34.9995838) ( 30d 0' 1.50"E, 34d59'58.50"N)
Center      (  10.0000005,  45.0000005) ( 10d 0' 0.00"E, 45d 0' 0.00"N)
Band 1 Block=128x128 Type=Int16, ColorInterp=Gray
  NoData Value=-32768

when running gdalinfo srtm.vft. Now, when I try to "query" the data, such as e.g. gdallocationinfo srtm.vft <your x coord> <your y coord>, with "appropriate" coordinates x and y, I get stuff like

Report:
  Location: (0P,18000L)
  Band 1:
    <LocationInfo><File>mosaic.tif</File></LocationInfo>
    Value: -32768

I tried different switches for gdallocationinfo, such as -wgs84 etc., to no avail.

My goal, as described in the related SE-question above, is to have a file with x,y,z, where x,y is a "locus" in the selected area, and z is the elevation of that locus. I also would like to sample the data such that any two pairs of coordinates specify locations that are "sufficiently" apart (say, 300-500 km).

The question, therefore, is: How exactly to query the above data (using gdallocationinfo or otherwise) to achieve my goal?

EDIT: I tried a different tool, gdal_translate: gdal_translate -of XYZ mosaik.tif geotiffas.xyz, after running gdalwarp *.tif mosaik.tif. But goetiffas.xyz contains nodata values only, i.e. -32768, as the z-coordinates.

2

Here's how it might work using only two gdal utilities:

Two initial comments.

  1. In my earlier suggestion (in the link above) I did not suggest to create a mosaic. You can avoid taking up that disk space by just using a vrt.
  2. I skipped one crucial parameter to the gdallocationinfo command: -geoloc. This tells gdal that the X and Y are in the data coordinate system (not row column numbers)

So I have a directory of several srtm (3 arcsec) data tiles, with names like n29_e034_3arc_v2.tif. Here are the commands and the output:

micha@sid-work:SRTM-3arc_v3$ gdalbuildvrt srtm.vrt n*.tif
0...10...20...30...40...50...60...70...80...90...100 - done.
micha@sid-work:SRTM-3arc_v3$ gdalinfo srtm.vrt 
Driver: VRT/Virtual Raster
Files: srtm.vrt
       n28_e032_3arc_v2.tif
       .... (many more rows) ....
       n33_e035_3arc_v2.tif
       n33_e036_3arc_v2.tif
Size is 6001, 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 = (31.999583333333334,34.000416666666666)
Pixel Size = (0.000833333333333,-0.000833333333333)
Corner Coordinates:
Upper Left  (  31.9995833,  34.0004167) ( 31d59'58.50"E, 34d 0' 1.50"N)
Lower Left  (  31.9995833,  27.9995833) ( 31d59'58.50"E, 27d59'58.50"N)
Upper Right (  37.0004167,  34.0004167) ( 37d 0' 1.50"E, 34d 0' 1.50"N)
Lower Right (  37.0004167,  27.9995833) ( 37d 0' 1.50"E, 27d59'58.50"N)
Center      (  34.5000000,  31.0000000) ( 34d30' 0.00"E, 31d 0' 0.00"N)
Band 1 Block=128x128 Type=Int16, ColorInterp=Gray
  NoData Value=-32767

The vrt covers from 32 deg (west side) to 37 deg, and 28 (south edge) to 34. Since you are interested in elevations at 500 km intervals, which is about 1/2 degree, you need to extract pixel values at 1/2 degree steps. Here's a loop (in a bash shell in linux) that can do that. Make sure to add the -geoloc option:

for x in 33 33.5 34 34.5 35 35.5; do for y in 28 28.5 29 29.5 30; do gdallocationinfo -geoloc srtm.vrt $x $y; done; done
.... (many more rows of output) ...
Report:
  Location: (4200P,6000L)
  Band 1:
    <LocationInfo><File>n28_e035_3arc_v2.tif</File><File>n29_e035_3arc_v2.tif</File></LocationInfo>
    Value: 1109
Report:
  Location: (4200P,5400L)
  Band 1:
    <LocationInfo><File>n29_e035_3arc_v2.tif</File></LocationInfo>
    Value: 1037
Report:
  Location: (4200P,4800L)
  Band 1:
    <LocationInfo><File>n29_e035_3arc_v2.tif</File><File>n30_e035_3arc_v2.tif</File></LocationInfo>
    Value: 1578

HTH

  • Fantastic, I'll give this a try now. Importantly, this avoids projecting the data to Cartesian coordinates. – Ilonpilaaja Apr 1 at 14:07

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.