Changeset 52063


Ignore:
Timestamp:
Jun 13, 2012 9:13:28 PM (4 years ago)
Author:
hamish
Message:

the guts of the tile snapping, download, and patching logic

File:
1 edited

Legend:

Unmodified
Added
Removed
  • grass-addons/grass6/raster/r.in.onearth/r.in.onearth.twms

    r52057 r52063  
    66# and Markus Neteler. Support for pre-tiled WMS server by Hamish Bowman     #
    77#                                                                           #
    8 # COPYRIGHT:    (C) 2005-2012 by the GRASS Development Team                 #
     8# COPYRIGHT:        (C) 2005-2012 by the GRASS Development Team             #
    99#                                                                           #
    10 #               This program is free software under the GNU General Public  #
    11 #               License (>=v2). Read the file COPYING that comes with GRASS #
    12 #               for details.                                                #
     10#            This program is free software under the GNU General Public     #
     11#            License (>=v2). Read the file COPYING that comes with GRASS    #
     12#            for details.                                                   #
    1313#                                                                           #
    1414#############################################################################
     
    140140
    141141eval `g.region -g`
     142
    142143if [ "$nsres" != "$ewres" ] ; then
    143144    g.message -e "East-West and North-South region resolution must be the same"
     
    147148fi
    148149
    149 case "$res" in
    150    0.125 | 0.25 | 0.5 | 1 | 2 | 4 | 8 | 16 | 32 | 64)
    151         value="valid"
    152         ;;
    153    *)
    154         value="invalid"
    155         g.message -e "Resolution must be one of "
    156         g.message -e "       {0.125 0.25 0.5 1 2 4 8 16 32 64}"
    157         exit 1
    158         ;;
    159 esac
    160 
    161 
    162 
    163 #######################################################
    164 # start of todo code dump
     150#### find appropriate tile resolution
     151fixed_res="
     1520.000244140625
     1530.00048828125
     1540.0009765625
     1550.001953125
     1560.00390625
     1570.0078125
     1580.015625
     1590.03125
     1600.0625
     1610.125
     1620.25
     163"
     164i=-3
     165TAKE=9999
     166TAKE_i=9999
     167MIN=9999
     168for val in $fixed_res ; do
     169   # find nearest resolution:
     170   #DIFF=`echo "$res $val" | awk '{printf("%.15g", $1 > $2 ? $1 - $2 : $2 - $1)}'`
     171   # find nearest finer resolution:
     172   DIFF=`echo "$res $val" | awk '{printf("%.15g", $1 > $2 ? $1 - $2 : 9999)}'`
     173   NEW_MIN_DIFF=`echo "$MIN $DIFF" | awk '{printf("%d", $1 < $2 ? 0 : 1)}'`
     174   if [ "$NEW_MIN_DIFF" -eq 1 ] ; then
     175      MIN="$DIFF"
     176      TAKE="$val"
     177      TAKE_i="$i"
     178   fi
     179   i=`expr $i + 1`
     180done
     181tile_res=`echo "$TAKE_i" | awk '{printf("%.15g", 2^$1)}'`
     182g.message -d message="min_diff=$MIN  [2^$TAKE_i = $tile_res  ($TAKE * 512)]"
     183
    165184
    166185
     
    172191
    173192
     193g.region res="$TAKE" -a
     194eval `g.region -g`
    174195
    175196# snap lat,long to upper-left corner of nearest tile grid node
    176197# Round lat to nearest grid node
    177 #todo: awkify:
    178 lat = 90.0 - round((90.0 - lat) / $res) * $res;
    179 lon = -180.0 + floor((180.0 + lon) / $res) * $res;
    180 
    181 #maybe floor() lower left lat/lon, ceil() for upper right lat/lon
    182 # ?
    183 
     198#lat = 90.0   - round( (90.0 - lat) / $res) * $res;
     199#lon = -180.0 + floor( (180.0 + lon) / $res) * $res;
     200
     201# diff is always positive here, so int() acts like floor().
     202n=`echo "$n $tile_res" | \
     203   awk '{ printf("%.15g", 90.0 - int((90.0 - $1) / $2) * $2) }'`
     204s=`echo "$s $tile_res" | awk '
     205 function ceil(x)
     206 {
     207   return (x == int(x)) ? x : int(x)+1
     208 }
     209 {
     210   printf("%.15g", 90.0 - ceil((90.0 - $1) / $2) * $2)
     211 }'`
     212
     213e=`echo "$e $tile_res" | \
     214   awk '{ printf("%.15g", -180.0 + int((180.0 + $1) / $2) * $2) }'`
     215w=`echo "$w $tile_res" | awk '
     216 function ceil(x)
     217 {
     218   return (x == int(x)) ? x : int(x)+1
     219 }
     220 {
     221   printf("%.15g", -180.0 + ceil((180.0 + $1) / $2) * $2)
     222 }'`
     223
     224#needed?
     225#g.region n=$n s=$s e=$e w=$w
     226
     227NUM_TILES_x=`echo "$e $w $tile_res" | awk '{ printf("%d", ($1 - $2) / $3) }'`
     228NUM_TILES_y=`echo "$n $s $tile_res" | awk '{ printf("%d", ($1 - $2) / $3) }'`
     229
     230TOTAL_TILES=`expr $NUM_TILES_x \* $NUM_TILES_y`
     231
     232if [ "$TOTAL_TILES" -gt 100 ] ; then
     233    g.message -w "You've told it to download $TOTAL_TILES tiles from NASA's server. Are you sure you want to do that?"
     234    read result
     235    case $result in
     236       y | yes | Y | Yes | YES)
     237          ;;
     238       *)
     239          g.message -i "Aborting."
     240          exit 1
     241          ;;
     242     esac
     243fi
     244
     245
     246#######################################################
     247# start of todo code dump
     248
     249
     250LAYER="global_mosaic_base"
     251STYLE="visual"
    184252
    185253BASE_URL="http://onearth.jpl.nasa.gov/wms.cgi"
    186254BASE_REQUEST="request=GetMap&layers=$LAYER&srs=EPSG:4326&width=512&height=512"
    187255BASE_META="format=image/jpeg&version=1.1.1&styles=$STYLE"
    188 URL="$BASE_URL?$BASE_REQUEST&$BBOX&$BASE_META"
    189 
    190 #bbox=minx,miny,maxx,maxy
    191 
    192 #top row
    193 minx=$(echo "$lon $res" | awk '{printf("%.9f", $1 - $2)}')
    194 miny="$lat"
    195 maxx="$lon"
    196 maxy=$(echo "$lat $res" | awk '{printf("%.9f", $1 + $2)}')
    197 BBOX_11="$minx,$miny,$maxx,$maxy"
    198 
    199 minx="$lon"
    200 miny="$lat"
    201 maxx=$(echo "$lon $res" | awk '{printf("%.9f", $1 + $2)}')
    202 maxy=$(echo "$lat $res" | awk '{printf("%.9f", $1 + $2)}')
    203 BBOX_21="$minx,$miny,$maxx,$maxy"
    204 
    205 minx=$(echo "$lon $res" | awk '{printf("%.9f", $1 + $2)}')
    206 miny="$lat"
    207 maxx=$(echo "$lon $res" | awk '{printf("%.9f", $1 + (2 * $2))}')
    208 maxy=$(echo "$lat $res" | awk '{printf("%.9f", $1 + $2)}')
    209 BBOX_31="$minx,$miny,$maxx,$maxy"
    210 
    211 #bottom row
    212 minx=$(echo "$lon $res" | awk '{printf("%.9f", $1 - $2)}')
    213 miny=$(echo "$lat $res" | awk '{printf("%.9f", $1 - $2)}')
    214 maxx="$lon"
    215 maxy="$lat"
    216 BBOX_12="$minx,$miny,$maxx,$maxy"
    217 
    218 minx="$lon"
    219 miny=$(echo "$lat $res" | awk '{printf("%.9f", $1 - $2)}')
    220 maxx=$(echo "$lon $res" | awk '{printf("%.9f", $1 + $2)}')
    221 maxy="$lat"
    222 BBOX_22="$minx,$miny,$maxx,$maxy"
    223 
    224 minx=$(echo "$lon $res" | awk '{printf("%.9f", $1 + $2)}')
    225 miny=$(echo "$lat $res" | awk '{printf("%.9f", $1 - $2)}')
    226 maxx=$(echo "$lon $res" | awk '{printf("%.9f", $1 + (2 * $2))}')
    227 maxy="$lat"
    228 BBOX_32="$minx,$miny,$maxx,$maxy"
    229 
    230 
    231 
    232 #top row
    233 URL="$BASE_URL?$BASE_REQUEST&bbox=$BBOX_11&$BASE_META"
    234 wget -nv "$URL" -O twms_11.jpg &
    235 
    236 URL="$BASE_URL?$BASE_REQUEST&bbox=$BBOX_21&$BASE_META"
    237 wget -nv "$URL" -O twms_21.jpg &
    238 
    239 URL="$BASE_URL?$BASE_REQUEST&bbox=$BBOX_31&$BASE_META"
    240 wget -nv "$URL" -O twms_31.jpg
    241 wait
    242 
    243 
    244 #bottom row
    245 URL="$BASE_URL?$BASE_REQUEST&bbox=$BBOX_12&$BASE_META"
    246 wget -nv "$URL" -O twms_12.jpg &
    247 
    248 URL="$BASE_URL?$BASE_REQUEST&bbox=$BBOX_22&$BASE_META"
    249 wget -nv "$URL" -O twms_22.jpg &
    250 
    251 URL="$BASE_URL?$BASE_REQUEST&bbox=$BBOX_32&$BASE_META"
    252 wget -nv "$URL" -O twms_32.jpg
    253 wait
     256
     257
     258# try a few at once, but don't get abusive
     259MAX_DL_JOBS=4
     260
     261for i in `seq $NUM_TILES_x` ; do
     262   i_str=`echo "$i" | awk '{printf("%04d", $1)}'`
     263   for j in `seq $NUM_TILES_y` ; do
     264      j_str=`echo "$j" | awk '{printf("%04d", $1)}'`
     265      minx=$(echo "$w $i $tile_res" | awk '{printf("%.9f", $1 + (($2 - 1) * $3))}')
     266      maxx=$(echo "$w $i $tile_res" | awk '{printf("%.9f", $1 + ($2 * $3))}')
     267      maxy=$(echo "$n $j $tile_res" | awk '{printf("%.9f", $1 - (($2 - 1) * $3))}')
     268      miny=$(echo "$n $j $tile_res" | awk '{printf("%.9f", $1 - ($2 * $3))}')
     269
     270      BBOX_TILE="$minx,$miny,$maxx,$maxy"
     271      BBOX="bbox=$BBOX_TILE"
     272      g.message -d message="  $BBOX"
     273
     274      URL="$BASE_URL?$BASE_REQUEST&$BBOX&$BASE_META"
     275
     276      CMD="wget -nv \"$URL\" -O twms_${i_str}_${j_str}.jpg"
     277
     278      MODULUS=`echo "$j $MAX_DL_JOBS" | awk '{print $1 % $2}'`
     279      g.message -d message="modulus=$MODULUS"
     280      if [ "$MODULUS" = "0" ] || [ "$j" = "$NUM_TILES_y" ] ; then
     281          # stall to let the background jobs finish
     282         g.message -v "Downloading tile [$i,$j] of [${NUM_TILES_x}x$NUM_TILES_y]."
     283         eval $CMD
     284         sleep 0.08
     285         wait
     286      else
     287         g.message -v "Downloading tile [$i,$j] of [${NUM_TILES_x}x$NUM_TILES_y],"
     288         eval $CMD &
     289         sleep 0.08
     290      fi
     291      wait
     292   done
     293done
    254294
    255295
    256296#### some checks before going on
    257 if [ `ls -1 twms_*.jpg | wc -l` -lt 6 ] ; then
    258     echo "ERROR: Tile(s) appear to be missing." 1>&2
    259     if [ "$verbose" = "true" ] ; then
    260         ls twms_*.jpg
    261     fi
    262     # todo: look to see which one is missing and try again.
    263     cd ..
    264     rm -rf "$TEMPDIR2"
     297if [ `ls -1 twms_*.jpg | wc -l` -lt "$TOTAL_TILES" ] ; then
     298    g.message -e "Tile(s) appear to be missing."
    265299    exit 1
    266300fi
     
    268302for file in twms_*.jpg ; do
    269303    if [ ! -s "$file" ] ; then
    270         echo "ERROR: <$file> appears to be empty." 1>&2
    271         if [ "$verbose" = "true" ] ; then
    272             ls -l twms_*.jpg
    273         fi
    274         cd ..
    275         rm -rf "$TEMPDIR2"
    276         exit 1
     304        echo "ERROR: <$file> appears to be empty." 1>&2
     305        exit 1
    277306    fi
    278307    if [ `file "$file" | grep -c JPEG` -eq 0 ] ; then
    279         echo "ERROR: <$file> appears to be bogus." 1>&2
    280         if [ `file "$file" | grep -c XML` -eq 1 ] && [ "$verbose" = "true" ] ; then
    281            cat "$file"
    282         fi
    283         cd ..
    284         rm -rf "$TEMPDIR2"
    285         exit 1
     308        echo "ERROR: <$file> appears to be bogus." 1>&2
     309        if [ `file "$file" | grep -c XML` -eq 1 ] && [ "$GRASS_VERBOSE" = "true" ] ; then
     310           cat "$file"
     311        fi
     312        exit 1
    286313    fi
    287314done
    288315
    289 if [ "$verbose" = "true" ] ; then
    290    echo "Converting to pnm ..." 1>&2
    291 fi
    292 
     316g.message -v "Converting to pnm ..." 1>&2
    293317for file in twms_*.jpg ; do
    294318    jpegtopnm "$file" > `basename "$file" .jpg`.pnm
     
    296320
    297321
    298 if [ "$verbose" = "true" ] ; then
    299    echo "Patching ..." 1>&2
    300 fi
    301 pnmcat -lr twms_11.pnm twms_21.pnm twms_31.pnm > row1.pnm
    302 pnmcat -lr twms_12.pnm twms_22.pnm twms_32.pnm > row2.pnm
    303 pnmcat -tb row1.pnm row2.pnm | pnmcut 0 0 1280 1024 | \
    304    pnmtojpeg --quality=75 > mosaic.jpg
    305 
    306 rm -f twms_*.pnm twms_*.jpg row[0-9].pnm
     322
     323g.message -v "Patching ..." 1>&2
     324
     325for j in `seq $NUM_TILES_y` ; do
     326   j_str=`echo "$j" | awk '{printf("%04d", $1)}'`
     327   tilefiles=`ls twms_*_${j_str}.pnm`
     328   #echo $tilefiles
     329   pnmcat -lr $tilefiles > "row_${j_str}.pnm"
     330done
     331
     332pnmcat -tb row_*.pnm | pnmtojpeg --quality=75 > mosaic_$$.jpg
     333
     334rm -f twms_*.pnm twms_*.jpg row*.pnm
     335
     336gdal_translate mosaic_$$.jpg mosaic_$$.tif \
     337  -a_srs EPSG:4326 -a_ullr $w $n $e $s 
     338  #-co COMPRESS=DEFLATE
     339
     340r.in.gdal in="mosaic_$$.tif" out="tmp_rio_$$" \
     341   title="NASA OnEarth WMS $LAYER ($STYLE)"
     342
     343for COLOR in red green blue ; do
     344   r.support "tmp_rio_$$.$COLOR" \
     345      history="Data downloaded from NASA's OnEarth WMS server with r.in.onearth"
     346done
     347
     348# free up the disk space ASAP
     349\rm "mosaic_$$.tif"
    307350
    308351
     
    352395{
    353396    if [ "$USE_GDALWARP" = "true" ] ; then
    354         g.message -v "**** CONVERT DATA  ****"
    355         #create the new imagename
     397        g.message -v "**** CONVERT DATA  ****"
     398        #create the new imagename
    356399        IMAGEFILE_GDALWARP="$TMPDIR/Image_${LAYER}_${STYLE}_${HEIGHT}_${WIDTH}_gdalwarp"
    357400
    358         #convert the data to the current location, create Erdas Imagine Images (HFA)
     401           #convert the data to the current location, create Erdas Imagine Images (HFA)
    359402        gdalwarp -s_srs "$SRC" -t_srs "`g.proj -wf`" -of HFA \
    360403            "$IMAGEFILE" "$IMAGEFILE_GDALWARP"
    361404        if [ $? -ne 0 ] ; then
    362           g.message -i '!-------- CAN NOT CONVERT DATA --------!'
    363           g.message -i '!------------ WILL BREAK --------------!'
     405          g.message -i '!-------- CAN NOT CONVERT DATA --------!'
     406          g.message -i '!------------ WILL BREAK --------------!'
    364407          exitprocedure
    365408        fi
    366         g.message -v "**** DATA CONVERTED ****"
     409        g.message -v "**** DATA CONVERTED ****"
    367410        #remove the old image and convert the name
    368411        rm -f "$IMAGEFILE"
     
    418461                g.message "Message from Server $NASASERVER"
    419462                echo " "
    420                 echo '!------------BEGIN-ERROR-MESSGAE--------------!'
     463                echo '!------------BEGIN-ERROR-MESSGAE--------------!'
    421464                cat "$IMAGEFILE"
    422                 echo '!-------------END-ERROR-MESSGAE---------------!'
     465                echo '!-------------END-ERROR-MESSGAE---------------!'
    423466            fi
    424467        fi
     
    451494if [ "$GIS_FLAG_F" -eq 1 ] ; then
    452495    if [ -z "$GIS_OPT_FILE"] ; then
    453         g.message -e "Please specify the output filename"
     496          g.message -e "Please specify the output filename"
    454497        exit 1
    455498    fi
     
    503546if [ $? -ne 0 ] && [ "$USE_GDALWARP" = "false" ] ; then
    504547    g.message -e "NASA OnEarth data are in Latitude/Longitude. The \
    505                   current location projection differs and you don't \
     548                  current location projection differs and you don't \
    506549                  have gdalwarp! STOP."
    507550    exit 1
Note: See TracChangeset for help on using the changeset viewer.