Raster Cookbook

From NGDCWiki

Revision as of 11:45, 8 February 2012 by Jcc (Talk | contribs)
(diff) ←Older revision | Current revision (diff) | Newer revision→ (diff)
Jump to: navigation, search

Contents

Cookbook for working with raster data in SDE

[Return to the Spatial Cookbook]

Use a IDRISI img file with ArcGIS

there are two important parts to this file, a .doc file and a .img. The former contains the metadata and the latter is a binary image file. Given two files, ekform.img and ekform.doc, an extract from the doc file looks like:

file title  : EPA Vegetation Form based on Kuchler
data type   : byte
file type   : binary
columns     : 940
rows        : 590
ref. system : albersus
ref. units  : m
unit dist.  : 1.0000000
min. X      : -2400000.0000000
max. X      : 2300000.0000000
min. Y      : 250000.0000000
max. Y      : 3200000.0000000
... snip ...

create a hdr file

a hdr file named in this case "ekform.hdr" is required. It's a simple parameter/value format with the values derived from the original doc file

ncols 940
nrows 590
nbands 1
nbits 8
byteorder I
layout bil
skipbytes 0
ulxmap -2400000.0000000
ulymap  3200000.0000000
xdim 1.0
ydim 1.0

in this case, the nbits (number of bits/cell) is 8 since the original file is declared as byte. The xdim, ydim correspond to the unit distance. The ulxmap, ulymap correspond to the "min. X", "max. Y" respectively

rename the img file

rename the original img file, in this case ekform.img becomes ekform.bil

ArcGIS should now be able to use the image if the .img and .hdr files are named with the same prefix and in the same directory

adjust metadata (optional)

once ArcGIS can read the image, you can choose to associate a projection, and apply null to background fill values using it.

Load a binary integer raster into SDE

  • Given an image file with the following specifications:

file name: lctmp12.img
cols: 720
rows: 360
data type: integer (16 bit)
file size = 360 * 720 * 2 = 518,400 bytes
xMin = -180.0
yMin = -90.0
xMax = 180.0
yMax = 90.0
resolution: 0.5 degrees
flag value: -999
min. value  : -533
max. value  : 417
assuming the xMin and yMin coordinates are for the lower left corner of the raster

  • create a header file in a text editor with the image name and a .hdr extension (i.e. "lctmp12.hdr")
ncols 720
nrows 360
nbands 1
nbits 16
byteorder I
layout bsq
skipbytes 0
ulxmap -180
ulymap  90
xdim 0.5
ydim 0.5

By adding the hdr file, the image can be displayed in ArcMap. However, because these are are signed integer data, they must be converted. (Images are assumed to have unsigned integers.)

  • Rename the data file to have a .bsq extension "lctmp12.bsq"

This is a single band image so it really does not matter whether it is named BIL or BSQ, but if the image is actually bil or bsq, name the file using the correct extension, and define the correct layout type in the header file. This step is not necessary if using the context menu in ArcCatalog

  • import the raster to an ESRI GRID
Arc: imagegrid lctmp12.bsq t bilinear
  • convert the data to unsigned integers (see <A HREF="#CASE9">Converting signed integer data to unsigned integer</A>)
Grid: t2 = con(t >= 32768, t - 65536, t)

values now range from -999 to 417

  • set the flag value to NoData
Grid: lctmp12 = setnull(t2 == -999, t2)

Now the Grid has the correct values and can be imported into SDE

Although sde_raster allows the NoData value to be specified during import, it does not support the Grid format.

  • import the Grid into SDE

use the context menu in ArcCatalog, "export -> raster to geodatabase". Note that this actually creates an "embedded raster catalog" with a single raster.

  • Grant privileges if necessary

context menu in ArcCatalog or sde_raster




Load an ascii float raster into SDE

  • Given a file with the following specifications:,br>

filename: area_co.txt
format: ASCII
datatype: float
One value per line
Raster size is 1332 rows by 1008 columns
The SW corner of the raster is -2736000,-2088000 meters (Lambert projection)
Cell size: 4km

  • file excerpt:
 1.404E+00
 9.368E-01
 6.705E-01
 7.266E-01
  • Insert header info into the start of the file:
ncols 1332
nrows 1008
xllcorner -2736000
yllcorner -2088000
cellsize 4000
  • Convert the text file into a ESRI Grid named "area_co"
Arc: asciigrid area_co.txt area_co float
  • (optional) Set flag values to NoData
  • (optional) Define the projection of the newly created grid
  • Import grid into SDE

use context menu in ArcCatalog, "export -> raster to geodatabase". Note that this actually creates a "embeded raster catalog" with a single raster.

  • Grant privileges if necessary

context menu in ArcCatalog or sde_raster




Add a color table to an SDE raster

  • Check for existing colortable

Check for its presence in ArcCatalog or by examining the SDE_AUX_xx table where "xx" indicates the SDE.RASTER_COLUMNS.RASTERCOLUMN_ID value.

SQL> select * from sde_aux_57;

RASTERBAND_ID       TYPE OBJECT
------------- ---------- ------
            1          2 0
            1          3 0

The record with a type "2" is the colortable record and the colortable itself stored in the LONG RAW column "OBJECT"

  • Drop the existing colortable with sderaster:
sderaster -o colormap -d -l wrzsoil,image -v 1 ^
-s cheetah-u [username] -p [password]

The argument "-v" refers to the "raster_id" and is usually "1". It is NOT the RASTERCOLUMN_ID value.

  • Add a new colortable

Use a TIFF image that contains the desired colortable. The contents of the image do not matter for this purpose, only the color table. So, one can use the same TIFF to provide an identical color table to multiple SDE rasters.

sderaster -o colormap -f wrzsoil.tif -l wrzsoil,image -v 1 ^ 
-s cheetah -u [username] -p [password]
  • Rebuild pyramids

Once the colortable is created, it is usually a good idea to rebuild the pyramids for the raster.

sderaster -o pyramid -l wrzsoil,image -v 1 -L -1 -I bilinear ^
-s cheetah -u [username] -p [password]
  • "-I bilinear" specifies bilinear resampling, "-L" indicates pyramid level and a value of "-1" creates a full pyramid (all levels)




Add a color table to a ESRI GRID

  • Externally create an CLR file for ESRI GRIDs.

This file has the format:

val0 R G B
val1 R G B
...

There must be a line in the file for each unique value in the GRID.
The file must be named <GRIDNAME>.clr and placed in the workspace directory.
The rrd and aux files may have to be removed and rebuilt for the new color table to be recognized.



preparing a colortable to add to a raster

Edit/Examine the colortable of a SDE raster

  • Export the raster into a BSQ file using sderaster
  • The ".clr" file contains the colortable. Edit with text editor
  • Upload the modified color table using sderaster



Storing the raster legend descriptions in a database

see ESRI Knowledgebase article 23728

  • Create a table
create table WRZSOIL_CLR (COLORMAP_INDEX VARCHAR2(3),DESCRIPTION VARCHAR2(64));
  • Populate the colortable and description values:
insert into WRZSOIL_CLR values ('0','Water/Ocean/Lake');
insert into WRZSOIL_CLR values ('1','Ferric Acrisol');
insert into WRZSOIL_CLR values ('2','Gleyic Acrisol');
  • Add/modify the layer in the AXL map configuration file
      <LAYER type="image" name="wrzsoil" visible="true" id="0">
         <DATASET name="GED.WRZSOIL.IMAGE" description="GED.WRZSOIL_CLR"
                  showcolormaplegend="true" workspace="sde_ws-0" />
         <IMAGEPROPERTIES  transcolor="0,0,0" />
      </LAYER>

Note that in the ArcXML 4.0.1 DTD, the "description" attribute for the <DATASET> element is undocumented and not listed . This means that it will not validate against the DTD unless you modify the DTD.



Copying a raster between SDE instances

One cannot use cut-and-paste in ArcCatalog to copy rasters (as you can vector layers). However, one can:

  • use ArcCatalog: Export -> Raster to Different Format -> Copy Raster
    • Specify configuration (DBTUNE) keyword, if desired
    • under Environment Settings/Raster Storage Settings, specify pyramid resampling technique, whether to calculate statistics, etc.
    • under Environment Settings/Raster Analysis Settings, specify cell size if desired
  • use ArcCatalog or sde_raster to create image file on disk and re-import into new instance
    • sde_raster can be used to adjust the pyramids after the copy if necessary.



Converting signed integer data to unsigned integer

  • Use the GRID module in ArcInfo Workstation or use ArcToolbox -> Spatial Analyst Tools -> Map Algebra:
newgrid = con(oldgrid >= 32768, oldgrid - 65536, oldgrid)



Set a flag value to NODATA

newgrid = setnull(oldgrid == -999, oldgrid)

loading a 16-bit signed integer raster (e.g. GLOBE data)

  • given a file "carib.bil" containing a single band of 16-bit signed integers, create a "carib.hdr" file like:
BYTEORDER       M
LAYOUT          BIL
NROWS           2160
NCOLS           2160
NBANDS          1
NBITS           16
BANDROWBYTES    4320
TOTALROWBYTES   4320
BANDGAPBYTES    0
NODATA          -500
ULXMAP          -85.99583330
ULYMAP          34.99583332
XDIM            0.00833333
YDIM            0.00833333
  • file size should be 2160 * 2160 * 2 = 9,331,200 bytes (rows * cols * bytes per cell)
  • use imagegrid to import into an ESRI GRID "carib":
imagegrid carib.bil carib
  • define the projection:
projectdefine grid carib
projection geographic
units dd
datum wgs84
zunits meters
parameters
  • adjust the grid to allow for signed integers:
carib2 = con(carib >= 32768, carib - 65536, carib)
  • set the nodata values. "-500" used as a flag value in GLOBE
carib3 = setnull(carib2 == -500, carib2)
  • use ArcCatalog to export to ArcSDE

load a TIFF into SDE using sderaster

  • Given a tiff image, wsiearth3.tif, w/colortable:
sderaster -o import -l wsiearth3,image -g -a 255,0,0 -f wsiearth3.tif ^
-s cheetah -u [username] -p [password]
  • the "-g" registers the image with the geodatabase. Without this, the image will not appear in the ArcCatalog/ArcMap table of contents
  • the "-a" indicates to convert all pixels with color 255,0,0 to NoData



Create a single GeoTIFF from the GLOBE tiles

  • download tiles and unzip them
  • rename tiles to have ".bil" extension
 for i in *.flt;do name=`echo $i |awk -F. '{print $1}'`; mv -i $name.flt $name.bil; done
  • gdal_merge.py -o globe.bil -of Ehdr -v *10g.bil
  • insure GDAL_DATA is set so that georef information available, e.g. export GDAL_DATA=$GDAL_HOME/data
  • gdal_translate -ot Int16 -of GTiff -a_srs EPSG:4326 -a_nodata -500 globe.bil globe.tif

Miscellenous Notes

  • creating statistics/analyzing a raster may change the appearance



References

Personal tools