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)]

edit ASCII file to have an azgcorr compatible header (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