Getting bands on Height/DEM/Raster? [closed]

Question

I can't seem to get rid of them, and I've no idea (OK scrap that, I have a vague idea) what's causing them.

How to reproduce

OK, so how did I end up with that.

• 1. I downloaded the ess-sh meridian 2 map from the OS Open data site
• 2. I unpacked the archive, and got the 'gridded height' folder

For the purposes of the image above, I only used the files in the NZ folder for now (I intend to do the entire UK once I perfect my method

• 3. I used the QGIS shape merge tool to combine all the shape files in the NZ folder into one large NZ grid of height points, each point is 200 meters apart (According to the measure tool)
• 4. I then loaded the combined shape file into QGIS, and did a visual inspection, each tile of points lined up perfectly, with no strange gaps (One tile starts 200 in from the edge of the tile, and the previous stops 200 after (That is it stops right on the grid square boundary (I've checked this using the OS grid locator shapes and everything is lined up on the boundaries using OSGB36)
• 5. I then used the QGIS interpolation tool, with the following settings:
• • Input layer : My NZ layer of gridded points
• • Interpolation attribute : Height
• • Interpolation method : TIN
• • Num of columns : 10000
• • Num of rows : 10000
• • Cell size X : 9.98000 (Auto)
• • Cell Size Y : 9.98000 (Auto)
• • X min, Y Min : 400100, 500100 (Min extents of entire NZ area)
• • X max, Y max : 499900, 599900 (Max extents of entire NZ area)
• 6. Hit OK and waited for the raster to render, the result you see above.

What I suspect

I suspect that the bands are some artefact of the shape-file joiner either not joining the source files correctly, or the data in the grids not being what I expect it to be.

However, IF there is discrepancies in how the OS has mapped these (and it wouldn't be the first time) then looking at the boundaries of the sources pre merge, I believe that I should also see these bands running vertically too.

---

I've since used a derivative of user30184's answer, which with a little help for here : http://darrencope.com/2010/05/07/merge-a-directory-of-shapefiles-using-ogr/ I turned into the following batch file:

@Echo off
Cls
if "%1" == "" goto NoParams

if not exist merged ( mkdir merged )

cd %1

for %%f in (*.shp) do (
if not exist ..\merged\%1merged.shp (
ogr2ogr -f "esri shapefile" ..\merged\%1merged.shp %%f -nln %1merged) else (
ogr2ogr -f "esri shapefile" -update -append ..\merged\%1merged.shp %%f -nln %1merged)
)

cd ..

goto Finished



:NoParams
ECHO This utility must be run from within the "Gridded_Height" folder as part of the OS     Meridian 2 map product
ECHO EG: merid2_essh_gb\data\Gridded_Height
ECHO.
ECHO When run, you must provide a 2 letter grid tile to process, the output of which     will be placed in a folder
ECHO called merged, EG: '%0 NZ'  will process all the shapes in the NZ tile square, and     produce
ECHO NZmerged.shp in the merged folder.
ECHO.

:Finished

Everything went smoothly, as anticipated, and so I ran the new shape file through the interpolation generator again, only to find that I still have these vertical lines on the raster.

As you can imagine, this is quite perplexing now, and I'm suspecting that either I have corrupted data (So am going to try and download a new copy) or ogr2ogr and other tools on the PC I'm using are broken (It is a possibility) so am going to try repeating the process on an entirely different machine and in a clean VM that's only ever had the geo tools I use installed.

For now however, still no dice, so all I can keep doing is fiddling, and see what happens.

I will keep this post updated in due course, if / when I get any further.

---

Further to the notes above, I've now taken the entire NZ layer and pushed the data into a PostGIS table, I've manually loaded each tile by hand, and then pulled this into QGIS and processed as described previously.

This has still produced the banding previously seen, but I now have some additional information.

Because I now have this in Postgres, I was able to overlay other layers onto the output in QGIS, doing this with the various grids I have from the Ordnance survey (I have a number of rectangular grids provided by them to aid positioning raster tiles), after trying a few different grids, Iv'e noted that these bands are exactly on the 10 km grid boundaries, but only vertically, horizontally works fine.

At this point I think I can rule out the shapefile merge tool in QGIS and any problems with Ogr2Ogr, this is now looking either like bad data (Which I have yet to find) or a problem in how the QGIS raster interpolation tool is producing the output.

Next step I think is going to be try the same procedure on a different machine, but accessing the same database data.

For reference (In case any one knows of any bugs) I'm using QGIS Valmiera (2.2.0)

in the process of loading all the shapes into individual tables inside Postgres, I came up with this script to merge them all into a single table (Posting it as it may be useful to others)

DO
$$
DECLARE
    tables        record;
    tablename     text;
    theschema     text;
    tile          text;
    tilenumber    text;
    sqlstring     text;
    
BEGIN

   SELECT 'temp' INTO theschema;

   FOR tables IN
      SELECT *
      FROM   pg_tables
      WHERE  schemaname = theschema
      ORDER  BY tablename
   LOOP
    
    SELECT quote_ident(tables.tablename) INTO tablename;
    SELECT substr(tables.tablename, 1 , 2) INTO tile;
    SELECT substr(tables.tablename, 3 , 2) INTO tilenumber;

    RAISE NOTICE 'Processing Table Name % (tile %, number %)', tablename, tile,     tilenumber;
    SELECT 'insert into meridian.gridded_height(originaltile, tilenumber, featurecode, height, geometry) (select ''' || tile || '''::character(2) as originaltile, ''' || tilenumber || '''::character(2) as tilenumber, "FEATURE_CO"::integer as featurecode, "HEIGHT" as height, geometry as geometry from ' || quote_ident(theschema) || '.' || tablename || ')' into sqlstring;
    EXECUTE sqlstring;
           
   END LOOP;
END$$;

The destination table schema the above script pushes the tile data into is:

CREATE TABLE meridian.gridded_height
(
  pkid serial NOT NULL,
  originaltile character(2),
  tilenumber character(2),
  featurecode integer,
  height double precision,
  geometry geometry(Point,27700),
  CONSTRAINT gridded_height_pkey PRIMARY KEY (pkid)
)
WITH (
  OIDS=FALSE

---

In follow up to user30184's most recent comment I'm using the following:

64 Bit Windows 7 with a 64 Bit build of QGIS, I don't have a separate install of GDAL (or any other geo utils, Ogr2Ogr as far as I'm aware runs out of my QGIS install)

Versions reported by QGIS (2.2.0 Valmiera) are as follows:

• QGIS version 2.2.0-Valmiera
• QGIS code revision c3a2817
• Compiled against Qt 4.8.5
• Running against Qt 4.8.5
• Compiled against GDAL/OGR 1.10.1
• Running against GDAL/OGR 1.10.1
• Compiled against GEOS 3.4.2-CAPI-1.8.2
• Running against GEOS 3.4.2-CAPI-1.8.2 r3921
• PostgreSQL Client Version 9.2.4
• SpatiaLite Version 4.1.1
• QWT Version 5.2.3
• PROJ.4 Version 480
• QScintilla2 Version 2.7.2

These versions have been taken directly from the QGIS about box

Further more, I've been doing some further analysis on the data.

The gridded height was apparently derived from the Landform panorama dataset sampled at 50 meters, however the distance between each point in the meridian grid appears to be 200 Meters.

When I overlay the 10 km grid (as previously mentioned) the grid lines run perfectly neatly between the grid points, looking at the co-ordinates of the height points, they are set in 100 meters from the edge of the grid, with 200 meter spacing.

This measures out at a 100 meter border between grid square edge and first column of points, but overall from the last column of the previous tile to the next, 200 meters to match the spacing for the rest correctly.

I've been trying different row and column counts in the interpolation generator, and have to some degree succeeded in removing some of the bands, and some combinations are worse than others, Iv'e not yet removed all of the bands, but given Iv'e found some changeable parameter that seems to be able to affect them, I believe that the data has not been produced in a fully suitable format by Ordnance survey to be used correctly for this purpose.

Following up on the new comments below, I'm going to try different versions of GDAL and co, see if that makes any difference at all.

Since last update, have also tried on a different machine, with same data and gotten the same results.

---

I've still not solved this problem, and I'm still no closer to why it occurred.

However, as time was running out for me and with help from a friend who knows ArcMap, I converted the points into greyscale DEMs using ArcMap with no problems at all.

I can't remember off the top of my head the exact sequence, but I do know it was much longer than the QGIS way of doing it.

I will when I get time however come back to this and continue to find a solution.

The following text was updated on the 10th July 2019 ==>

This question was originally asked when QGIS was in the version 2 branch, specifically version 2.2.

In recent years as QGIS has moved to the V3.x branch, I've not seen, this problem re-occur. This leads me to believe that it may have been caused by a bug in QGis 2.2

This issue as far as I've seen cannot be reproduced unless you are still using QGIS V2.2.

Answer 1

I can't say exactly why you got those lines but I can show how you can make build a good DEM. The weak point is probably the "Merge shapefiles to one" tool. Don't use it, use GDAL instead. The merge tool seems to force reprojecting into EPSG:4326 which is at least not desired always.

This is what I did

Take the first shapefile and write it into a new shapefile "total.shp"

Then take the second one and append it to "total.shp". Do the same for all the shapefiles. For the whole UK you'll want to do a script.

ogr2ogr -f "ESRI Shapefile" total.shp NZ00.SHP -nln total
ogr2ogr -f "ESRI Shapefile" -update -append total.shp NZ01.SHP -nln total
ogr2ogr -f "ESRI Shapefile" -update -append total.shp NZ02.SHP -nln total
ogr2ogr -f "ESRI Shapefile" -update -append total.shp NZ03.SHP -nln total
ogr2ogr -f "ESRI Shapefile" -update -append total.shp NZ04.SHP -nln total
ogr2ogr -f "ESRI Shapefile" -update -append total.shp NZ05.SHP -nln total
...

Make DEM as you did before with QGIS. At least for me the result looks good.