landfire.sh

   1 #!/bin/bash
   2
   3 output_dir="/home/jbeezley/wrfdata/geog/landfire_data"
   4
   5 function get_param
   6 {
   7   local fname
   8   local pname
   9   fname=$1
  10   pname=$2
  11   pval=$(grep -a $pname $fname | sed "s/^${pname}  *//" | sed 's/ *//g' | tr -d '\r' )
  12 }
  13
  14 function split_string
  15 {
  16   local str
  17   str=$(echo $1 | sed 's/^ *//' | sed 's/ *$//' | sed 's/  */ /g')
  18   pdeg=$(echo $str | cut -s -d " " -f 1)
  19   pmin=$(echo $str | cut -s -d " " -f 2)
  20   psec=$(echo $str | cut -s -d " " -f 3)
  21 }
  22
  23 function get_origin_params
  24 {
  25   local fname
  26   local pval
  27   local sstr
  28   fname=$1
  29   sstr='\/\* 1st standard parallel'
  30   par1str=$(grep -a "$sstr" $fname | sed "s/ *${sstr}.*$//" | tr -d '\r')
  31   split_string "$par1str"
  32   par1deg=$pdeg
  33   par1min=$pmin
  34   par1sec=$psec
  35   sstr='\/\* 2nd standard parallel'
  36   par2str=$(grep -a "$sstr" $fname | sed "s/ *${sstr}.*$//" | tr -d '\r')
  37   split_string "$par2str"
  38   par2deg=$pdeg
  39   par2min=$pmin
  40   par2sec=$psec
  41   sstr='\/\* central meridian'
  42   merstr=$(grep -a "$sstr" $fname | sed "s/ *${sstr}.*$//" | tr -d '\r')
  43   split_string "$merstr"
  44   merdeg=$pdeg
  45   mermin=$pmin
  46   mersec=$psec
  47   sstr='\/\* latitude of projection'
  48   latstr=$(grep -a "$sstr" $fname | sed "s/ *${sstr}.*$//" | tr -d '\r')
  49   split_string "$latstr"
  50   latdeg=$pdeg
  51   latmin=$pmin
  52   latsec=$psec
  53   sstr='\/\* false easting'
  54   falseeast=$(grep -a "$sstr" $fname | sed "s/ *${sstr}.*$//" | tr -d '\r')
  55   sstr='\/\* false northing'
  56   falsenorth=$(grep -a "$sstr" $fname | sed "s/ *${sstr}.*$//" | tr -d '\r')
  57 }
  58
  59 function read_headers
  60 {
  61   local dir fname
  62   dir=$1
  63   fname="${dir}/${dir}.prj"
  64   get_param $fname Projection
  65   proj="$pval"
  66   get_param $fname Datum
  67   datum="$pval"
  68   get_param $fname Spheroid
  69   speroid="$pval"
  70   get_param $fname Units
  71   units="$pval"
  72   get_param $fname Zunits
  73   zunits="$pval"
  74   get_param $fname Xshift
  75   xshift="$pval"
  76   get_param $fname Yshift
  77   yshift="$pval"
  78
  79   get_origin_params $fname
  80
  81   fname="${dir}/${dir}.hdr"
  82   get_param $fname ncols
  83   ncols="$pval"
  84   get_param $fname nrows
  85   nrows="$pval"
  86   get_param $fname xllcorner
  87   xll="$pval"
  88   get_param $fname yllcorner
  89   yll="$pval"
  90   get_param $fname cellsize
  91   cellsize="$pval"
  92   get_param $fname NODATA_value
  93   nodata="$pval"
  94   get_param $fname byteorder
  95   byteorder="$pval"
  96
  97 }
  98
  99 read_headers $1
 100 rm -fr "$output_dir" &> /dev/null
 101 mkdir -p "$output_dir" &> /dev/null
 102 idx="$output_dir/index"
 103
 104 if [ "$proj" == "ALBERS" ] && [ "$datum" == "NAD83" ] ; then
 105   nproj=1
 106 #  echo "projection=albers_nad83" >> $idx
 107   echo "projection=regular_ll" >> $idx
 108 else
 109   echo "Projection $proj $datum not supported"
 110   exit 1
 111 fi
 112
 113 echo 'type=categorical' >> $idx
 114 echo 'units="category"' >> $idx
 115 echo 'description="Anderson 13 Fire Behavior Fuel Models"' >> $idx
 116
 117 if [ $nproj -eq 1 ] && [ $units == METERS ] ; then
 118   echo "dx=$cellsize" >> $idx
 119   echo "dy=$cellsize" >> $idx
 120 else
 121   echo "Unit conversion for $units to meters not supported"
 122   exit 1
 123 fi
 124
 125 echo 'known_x=1' >> $idx
 126 echo 'known_y=1' >> $idx
 127
 128 echo "WARNING: reference long/lat is not implemented yet, setting to static value"
 129 echo 'known_lat=35.7337' >> $idx
 130 echo 'known_lon=-107.3452' >> $idx
 131
 132 echo 'wordsize=2' >> $idx
 133 echo 'category_min=1' >> $idx
 134 echo 'category_max=14' >> $idx
 135 echo 'tile_bdr=0' >> $idx
 136 echo 'missing_value=0' >> $idx
 137 echo 'scale_factor=1' >> $idx
 138 echo 'row_order=bottom_top' >> $idx
 139
 140 if [ $byteorder == LSBFIRST ] ; then
 141   echo 'endian=little' >> $idx
 142 else
 143   echo 'endian=big' >> $idx
 144 fi
 145
 146 echo "tile_x=$nrows" >> $idx
 147 echo "tile_y=$ncols" >> $idx
 148 echo "tile_z=14" >> $idx
 149
 150 ./convert_landfire.x ${1}/${1}.flt $ncols $nrows
 151 nfname="$output_dir/$(printf "%05d,%05d.%05d,%05d" 1 $ncols 1 $nrows)"
 152 mv "${1}/${1}.flt_tmp" "$nfname"
 153