Creation of DEMs from NextMap data

(fast notes)

load arc files into grass

• import->raster->gdal
• pick file
• choose create new location if this is the first import
• run (will fail for UK as it needs user input)
• copy and paste command into grass shell so you can respond (e.g. r.in.gdal -e input=/data/aegean1/nextmap_dems/neodc/by_tile/by_product/dsm/nk/nk02/nk02dsm/w001001.adf output=nk02
• if you need to import multiple rasters, you can use the location created, but must select extend region
• (need to figure out how to join rasters into one, then subset in grass)

export as ASCII grid

[you can load this grid into ENVI to test against vector layers, etc - file->open external->generic->ascii, then georeference (edit headers->mapinfo, etc)]

The azgcorr format is plain ASCII text, with decimal numbers separated by spaces, one row per row of the DEM. At the top of the file, you need a single line that describes the geographic position of the DEM.

Run the script to convert the ASCII file to have an azgcorr compatible header. Note this will delete the original ASCII file unless you specify -nodel at the end

# Run conversion script, input file is the first argument, output file is the second
demheaderconvert.sh asciifile.dem azgcorrfile.dem [-nodel]

If for any reason the script doesn't work, edit the ASCII file manually (first line must be as follows, then DEM data):

 -ed or c r xm ym xx yx gi : is the basic flat file ascii DEM definition
            : or = data order =0 rows S->N, =1 rows N->S
            : c = cols, r = rows, xm, ym = SW, xx, yx = NE corners
            : gi = grid increments; grid values separated by spaces
            : ** NB: for this definition georeferencing of the DEM grid is assumed
            :    to be at the CENTRE of the cell defined by the grid coordinates
            : ** if this is NOTthe case use option -edx, below

 e.g.
    1 2000 2000 400000 850000 410000 860000 5
 N->s Xsize YSize Xmin Ymin   Xmax   Ymax   grid size (5m for nextmap)

Use DEM with azgcorr -eh demfilename

---

Command line method

Starting from a valid location (in correct projection?)

 # read in a tile
r.in.gdal -e input=/data/aegean1/nextmap_dems/neodc/by_tile/by_product/dsm/tl/tl17/tl17dsm/w001001.adf output=tl17
 # set the region of interest to the (two) tiles we want to export
g.region rast=tl17,tl27
 # patch them together into a single raster
r.patch input=tl17,tl27 output=tl17and27
 # write out as an ASCII DEM (needs editing as above)
r.out.ascii input=tl17and27 output=/tmp/tl17and27.dem null=0

---

To add a zero-level area to a map (ie sea level):

In GRASS:

# Restrict active region to area you're interested in - must use eastings/northings, lat/long seems to have problems determining valid values
g.region n=601000 e=451000 s=579000 w=389000

# Check region is restricted correctly
g.region -gb

# Create vector map of active region, output param is name of output map
v.in.region output=v_sealevel

# Convert vector map to raster map
v.to.rast input=v_sealevel output=r_sealevel use=val value=0.0

# Create raster maps for areas with OS data as first line only in section above
r.in.gdal -e input=/data/aegean1/nextmap_dems/neodc/by_tile/by_product/dsm/ny/ny98/ny98dsm/w001001.adf output=ny98
r.in.gdal -e input=/data/aegean1/nextmap_dems/neodc/by_tile/by_product/dsm/ny/ny99/ny99dsm/w001001.adf output=ny99
...

# Patch raster maps together as section above
# NOTE SEA LEVEL MAP MUST BE LAST IN LIST, since Grass starts with a blank map and adds patched raster maps in order
r.patch input=ny98,ny99,...,r_sealevel output=map_with_added_sea_level

# Create DEM from raster map as section above
r.out.ascii input=map_with_added_sea_level output=/tmp/map_with_added_sea_level.dem null=0