Changeset 52185


Ignore:
Timestamp:
Jun 22, 2012 3:47:47 AM (4 years ago)
Author:
hamish
Message:

it's alive

File:
1 edited

Legend:

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

    r52176 r52185  
    22#############################################################################
    33# Download and import satellite images direct from the                      #
    4 # NASA onearth WMS server into GRASS.                                       #
    5 # written by Soeren Gebbert 11/2005 soerengebbert AT gmx de                 #
    6 # and Markus Neteler. Support for pre-tiled WMS server by Hamish Bowman     #
     4#  NASA onearth WMS server into GRASS.                                      #
     5# Written by Soeren Gebbert 11/2005 soerengebbert AT gmx de                 #
     6#  and Markus Neteler.                                                      #
     7# Rewritten to support for pre-tiled WMS server by Hamish Bowman            #
    78#                                                                           #
    89# COPYRIGHT:        (C) 2005-2012 by the GRASS Development Team             #
     
    1516
    1617#%Module
    17 #%  description: Download and import satellite images direct from the NASA onearth WMS server into GRASS or to a geo-tiff image file.
     18#%  description: Download and import satellite images direct from the NASA OnEarth Tiled WMS server into GRASS or to a GeoTIFF image file.
    1819#%End
    1920#%option
     
    2627#%option
    2728#%  key: file
    28 #%  gisprompt: file,file,file
     29#%  gisprompt: new_file,file,file
    2930#%  type: string
    3031#%  description: Output file name prefix
     
    139140
    140141
    141 ############
    142 # for now just support lat/lon
    143 USE_GDALWARP="false"
    144 ############
    145 
    146 
    147 
    148142# Get the data from the NASA server
    149143GetData()
    150144{
    151145   IMPORT="true" #default
     146
     147   # Compose the filename for the assembled GeoTIFF
     148   IMAGEFILE="$TMPDIR/OnEarth_${LAYER}_${STYLE}_$$.tif"
     149
     150   g.message -v "Download Data ..."
     151   g.message -v "Requesting Data from '$BASE_URL'"
     152
    152153   #local STRING="request=GetMap&layers=$LAYER&srs=$SRS&width=$WIDTH&height=$HEIGHT&bbox=$w,$s,$e,$n&format=$FORMAT&version=1.1.0&styles=$STYLE$TIME"
    153154   BASE_REQUEST="request=GetMap&layers=$LAYER&srs=$SRS&width=$TILE_SIZE&height=$TILE_SIZE"
    154155   BASE_META="format=$FORMAT&version=1.1.1&styles=$STYLE$TIME"
    155 
    156 
    157    local REQ_STRING="$BASE_REQUEST&$BBOX&$BASE_META"
    158    URL="$BASE_URL?$REQ_STRING"
    159 
    160    # Compose the filename
    161    IMAGEFILE="$TMPDIR/OnEarth_${LAYER}_${STYLE}_$$.tif"
    162    g.message -v "Download Data ..."
    163    g.message -v "Requesting Data from '$BASE_URL'"
     156   # what the request would look like if we had a general-purpose WMS server:
     157   #BBOX="bbox=$w,$s,$e,$n"
     158   #REQ_STRING="$BASE_REQUEST&$BBOX&$BASE_META"
     159   #WMS_URL="$BASE_URL?$REQ_STRING"
    164160
    165161   #download the File from the Server
    166162   #wget $WGET_OPTIONS --post-data="$REQ_STRING" "$BASE_URL" -O "$IMAGEFILE"
    167163
    168    OLDDIR=`pwd`
    169    cd "$TMDDIR"
    170    
    171164   for i in `seq $NUM_TILES_x` ; do
    172165      i_str=`echo "$i" | awk '{printf("%04d", $1)}'`
     166
     167      # avoid tall+thin jobs from running away
     168      MODULUS=`echo "$i $MAX_DL_JOBS" | awk '{print $1 % $2}'`
     169      if [ "$MODULUS" = "0" ] ; then
     170          wait
     171      fi
     172
    173173      for j in `seq $NUM_TILES_y` ; do
    174174         j_str=`echo "$j" | awk '{printf("%04d", $1)}'`
     
    182182         g.message -d message="  $BBOX"
    183183   
    184          local REQ_STRING="$BASE_REQUEST&$BBOX&$BASE_META"
    185          g.message -d "URL is [$BASE_URL?$REQ_STRING]"
     184         REQ_STRING="$BASE_REQUEST&$BBOX&$BASE_META"
     185         g.message -d message="URL is [$BASE_URL?$REQ_STRING]"
    186186   
    187187         #CMD="wget -nv \"$URL\" -O twms_${i_str}_${j_str}.jpg"
    188          CMD="wget -nv --post-data="$REQ_STRING" "$BASE_URL" -O twms_${i_str}_${j_str}.jpg"
    189    
    190    
     188         CMD="wget -nv --post-data='$REQ_STRING' '$BASE_URL' -O 'twms_${i_str}_${j_str}.jpg'"
     189
    191190         MODULUS=`echo "$j $MAX_DL_JOBS" | awk '{print $1 % $2}'`
    192191         g.message -d message="modulus=$MODULUS"
    193          if [ "$MODULUS" = "0" ] || [ "$j" = "$NUM_TILES_y" ] ; then
     192         if [ "$MODULUS" = "0" ] ; then
    194193             # stall to let the background jobs finish
    195             g.message -v "Downloading tile [$i,$j] of [${NUM_TILES_x}x$NUM_TILES_y]."
     194            g.message "Downloading tile [$i,$j] of [${NUM_TILES_x}x$NUM_TILES_y]."
    196195            eval $CMD
    197             EXITCODE=$?
    198             wait
    199             if [ "$EXITCODE" -ne 0 ] ; then
    200                g.message -e "wget was not able to download the data"
    201                IMPORT="false"
    202                return 1
     196            EXITCODE=$?
     197            wait
     198            if [ "$EXITCODE" -ne 0 ] ; then
     199               #retry
     200               sleep 1
     201               eval $CMD
     202               if [ "$EXITCODE" -ne 0 ] ; then
     203                  g.message -e "wget was not able to download the data"
     204                  IMPORT="false"
     205                  return 1
     206               fi
    203207            fi
    204208         else
    205             g.message -v "Downloading tile [$i,$j] of [${NUM_TILES_x}x$NUM_TILES_y],"
     209            g.message "Downloading tile [$i,$j] of [${NUM_TILES_x}x$NUM_TILES_y],"
    206210            eval $CMD &
    207211         fi
    208          wait
    209212      done
    210213   done
    211    
     214   wait
     215
    212216   #### some checks before going on
    213217   if [ `ls -1 twms_*.jpg | wc -l` -lt "$TOTAL_TILES" ] ; then
     
    218222   
    219223   for file in twms_*.jpg ; do
    220        if [ ! -s "$file" ] ; then
    221         g.message -e "<$file> appears to be empty."
    222         IMPORT="false"
    223         return 1
    224        fi
    225        if [ `file "$file" | grep -c JPEG` -eq 0 ] ; then
    226         g.message -e "<$file> appears to be bogus."
    227         if [ `file "$file" | grep -c XML` -eq 1 ] && [ "$GRASS_VERBOSE" = "true" ] ; then
    228            cat "$file"
    229         fi
    230         IMPORT="false"
    231         return 1
    232        fi
     224        if [ ! -s "$file" ] ; then
     225            g.message -e "<$file> appears to be empty."
     226            IMPORT="false"
     227            return 1
     228        fi
     229        if [ `file "$file" | grep -c JPEG` -eq 0 ] ; then
     230            g.message -e "<$file> appears to be bogus."
     231            if [ `file "$file" | grep -c XML` -eq 1 ] && [ "$GRASS_VERBOSE" = "true" ] ; then
     232                g.message message="--------------------------------------"
     233                cat "$file"
     234                g.message message="--------------------------------------"
     235            fi
     236            IMPORT="false"
     237            return 1
     238        fi
    233239   done
    234    
    235    g.message -v "Converting to pnm ..." 1>&2
     240
     241
     242   g.message -v "Converting to pnm ..."
    236243   for file in twms_*.jpg ; do
    237244       jpegtopnm "$file" > `basename "$file" .jpg`.pnm
    238245   done
    239246   
    240    g.message -v "Patching ..." 1>&2
     247   g.message -v "Patching ..."
    241248   
    242249   for j in `seq $NUM_TILES_y` ; do
     
    256263
    257264   # Convert to GeoTIFF
     265   # gdalwarp can get confused if the file already exists (or at least it used to)
     266   if [ -e "mosaic_$$.tif" ] ; then
     267       g.message -e "A file called [$TMPDIR/mosaic_$$.tif] already exists. Aborting."
     268       exit 1
     269   fi
    258270   gdal_translate "mosaic_$$.jpg" "mosaic_$$.tif" \
    259271      -a_srs EPSG:4326 -a_ullr $w $n $e $s
     
    280292{
    281293    if [ "$USE_GDALWARP" = "true" ] ; then
    282         g.message -v "**** CONVERT DATA  ****"
     294        g.message -v "Reprojecting the image ..."
    283295        #create the new imagename
    284         IMAGEFILE_GDALWARP="$TMPDIR/Image_${LAYER}_${STYLE}_${HEIGHT}_${WIDTH}_gdalwarp"
    285 
    286            #convert the data to the current location, create Erdas Imagine Images (HFA)
    287         gdalwarp -s_srs "$SRS" -t_srs "`g.proj -wf`" -of HFA \
    288             "$IMAGEFILE" "$IMAGEFILE_GDALWARP"
     296        IMAGE_WARPED="$TMPDIR/`basename "$IMAGEFILE" .tif`_warped.tif"
     297
     298        #convert the data to the current location, create Erdas Imagine Images (HFA)
     299        gdalwarp -s_srs "$SRS" -t_srs "`g.proj -wf`" \
     300            -dstalpha "$IMAGEFILE" "$IMAGE_WARPED"
    289301        if [ $? -ne 0 ] ; then
    290           g.message -i '!-------- CAN NOT CONVERT DATA --------!'
    291           g.message -i '!------------ WILL BREAK --------------!'
     302          g.message -i 'Reprojection failed. Aborting.'
    292303          exitprocedure
    293304        fi
    294         g.message -v "**** DATA CONVERTED ****"
     305        g.message -v "Reprojection successful"
    295306        #remove the old image and convert the name
    296307        rm -f "$IMAGEFILE"
    297         IMAGEFILE="$IMAGEFILE_GDALWARP"
     308        IMAGEFILE="$IMAGE_WARPED"
    298309        return 0
    299310    fi
     
    321332
    322333    if [ $? -ne 0 ] ; then
    323         g.message -e "Downloaded file is not supported by gdal, or cannot be imported"
    324         #if [ $ReturnValueGdalBug -eq 0 ] ; then
    325         #    echo "There was a problem while downloading the file, maybe you should try it again."
    326         #fi
    327         #If the File is XML, then cat the contents to stdout
    328         echo "$FILETYPE" | grep XML > /dev/null
    329         if [ $? -eq 0 ] ; then
    330             g.message " "
    331             g.message "Message from Server $BASE_URL"
    332             echo " " 1>&2
    333             echo '!------------BEGIN-ERROR-MESSGAE--------------!' 1>&2
    334             cat "$IMAGEFILE" 1>&2
    335             echo '!-------------END-ERROR-MESSGAE---------------!' 1>&2
    336         fi
     334        g.message -e "Downloaded file is not supported by GDAL, or cannot be imported"
    337335    fi
    338336
    339337    g.message -d "Data checked & was ok."
     338
    340339
    341340    ### Copy or import
    342341    if [ "$GIS_FLAG_F" -eq 1 ] ; then
    343342        # Copy the data to the output file
    344         g.message -v "Creating output file $GIS_OPT_FILE$TYPE$STYLE$FILE_EXTENT"
    345         cp "$IMAGEFILE" "$GIS_OPT_FILE$TYPE$STYLE$FILE_EXTENT"
     343        OUTNAME=`echo "$GIS_OPT_OUTPUT" | sed -e 's/_$//'`
     344        OUTFILE="${OUTNAME}_${TYPE}_$STYLE.tif"
     345
     346        g.message -v "Creating output file [$OUTFILE]"
     347        cp -f "$IMAGEFILE" "$WORKDIR/$OUTFILE"
    346348    else
    347349        # Warp the data!
     
    353355        OUTNAME="${OUTNAME}_${TYPE}_$STYLE"
    354356
    355         r.in.gdal input="$IMAGEFILE" output="$OUTNAME" \
    356            title="NASA OnEarth WMS $LAYER ($STYLE)"
    357 
     357        r.in.gdal input="$IMAGEFILE" output="$OUTNAME.tmp_$$" --quiet
     358
     359        # crop away no data introduced by differing convergence angle
     360        g.region "rast=$OUTNAME.tmp_$$.alpha"
    358361        for COLOR in red green blue ; do
    359            r.support "$OUTNAME.$COLOR" \
    360               history="Data downloaded from NASA's OnEarth WMS server with r.in.onearth"
     362            r.mapcalc "$OUTNAME.$COLOR = \
     363              if($OUTNAME.tmp_$$.alpha == 0, null(), $OUTNAME.tmp_$$.$COLOR)" &
    361364        done
     365        wait
     366
     367        for COLOR in red green blue ; do
     368            r.colors "$OUTNAME.$COLOR" color=grey255 --quiet
     369            r.support "$OUTNAME.$COLOR" \
     370              title="NASA OnEarth WMS $LAYER ($STYLE)"
     371              source1="Data downloaded from NASA's OnEarth WMS server" \
     372              source2="$LAYER ($STYLE)" \
     373              description="generated by r.in.onearth"
     374        done
     375       
     376        g.mremove -f rast="$OUTNAME.tmp_$$.*" --quiet
    362377    fi
    363378
     
    424439
    425440
    426 
    427 #######################################################
    428 # start of old code
    429 
    430441#Now get the LatLong Boundingbox
    431 if [ "$IS_4326" = "false" ] && [ "$USE_GDALWARP" = "true" ] ; then
     442if [ "$IS_4326" = "true" ] ; then
     443    #We have LatLong projection, no warp is needed!
     444    USE_GDALWARP="false"
     445else
    432446    eval `g.region -gb`
    433447    n="$ll_n"
     
    435449    e="$ll_e"
    436450    w="$ll_w"
    437     g.message -v "LatLong wgs84 bounding box = N $n S $s W $w E $e"
    438 else
    439     #We have LatLong projection, no warp is needed!
    440     USE_GDALWARP="false"
    441 fi
    442 
    443 if [ 0 -eq 1 ] ; then
    444     #There is a bug in nasa WMS service, it provides images which are lager then
    445     #the world :(, we have to crop the images
    446     if [ "$n" = "90" -a "$s" = "-90" ] && \
    447        [ "$w" = "-180" -a "$e" = "180" ] ; then
    448 
    449         # check if we have bc
    450         if [ ! -x "`which bc`" ] ; then
    451             g.message -e "bc required, please install first"
    452             exit 1
    453         fi
    454         #We request a smaller image from the wms server
    455         n=`echo "$n - 0.001" | bc`
    456         s=`echo "$s + 0.001" | bc`
    457         e=`echo "$e - 0.001" | bc`
    458         w=`echo "$w + 0.001" | bc`
    459     fi
    460 fi
    461 ###old
    462 
    463 # end of old code
    464 #######################################################
    465 
    466 
     451    g.message -d "Lat-Long WGS84 bounding box was = N $n S $s W $w E $e"
     452fi
    467453
    468454
     
    475461fi
    476462
    477 
    478 
    479 
    480 
    481 #######################################################
    482 # start of todo code dump
    483 
    484 eval `g.region -g`
     463eval `g.region -g | grep 'res='`
    485464
    486465if [ "$nsres" != "$ewres" ] ; then
     
    490469    res="$nsres"
    491470fi
     471
     472if [ $PROJ_CODE -ne 3 ] ; then
     473   # if not lat/lon, convert map units resolution to degrees
     474   TO_METERS=`g.proj -p | grep '^meters.*:' | awk '{print $3}'`
     475   res=`echo "$nsres $TO_METERS" | awk '{ printf("%.11f", $1 / ($2 * 1852.0 * 60))}'`
     476
     477fi
     478
    492479
    493480#### find appropriate tile resolution
     
    528515fi
    529516tile_res=`echo "$TAKE_i" | awk '{printf("%.15g", 2^$1)}'`
    530 g.message -i message="min_diff=$MIN  [2^$TAKE_i = $tile_res  ($TAKE * 512)]"
    531 
    532 
    533 #### adjust region to snap to grid
     517g.message -d message="min_diff=$MIN  [2^$TAKE_i = $tile_res  ($TAKE * 512)]"
     518
     519
     520#### adjust region to snap to fixed list of available resolutions
    534521# setup internal region
    535522g.region save="tmp_tiled_wms.$$"
     
    538525
    539526
    540 g.region res="$TAKE" -a
    541 eval `g.region -g`
    542 
     527### now snap the current region to the best-fit tiled grid
     528#d.region.box   #orig
     529if [ $PROJ_CODE -ne 3 ] ; then
     530   # if not lat/lon, convert map units resolution to degrees
     531   local_res=`echo "$TAKE $TO_METERS" | awk '{ printf("%.15g\n", $1 * $2 * 1852.0 * 60)}'`
     532   g.region res="$local_res" -a
     533   eval `g.region -gb`
     534   n="$ll_n"
     535   s="$ll_s"
     536   e="$ll_e"
     537   w="$ll_w"
     538   g.message -d "Lat-Long WGS84 bounding box aligned to = N $n S $s W $w E $e"
     539else
     540   g.region res="$TAKE" -a
     541   eval `g.region -g`
     542fi
     543#d.redraw
     544#d.region.box   #after resolution snap
     545
     546
     547
     548#### begin snap to nearest tile
    543549# snap lat,long to upper-left corner of nearest tile grid node
    544550# Round lat to nearest grid node
     
    569575 }'`
    570576
    571 # needed?
     577g.message -d "Lat-Long WGS84 bounding box snapped to = N $n S $s W $w E $e"
     578
     579# not needed:
    572580#g.region n=$n s=$s e=$e w=$w
     581#d.redraw
     582#d.region.box
     583
    573584
    574585#NUM_TILES_x=`echo "$e $w $tile_res" | awk '{ printf("%d", 1 + (($1 - $2) / $3)) }'`
     
    577588NUM_TILES_y=`echo "$n $s $tile_res" | awk '{ printf("%d", ($1 - $2) / $3) }'`
    578589TOTAL_TILES=`expr $NUM_TILES_x \* $NUM_TILES_y`
     590
     591g.message -v "There are $NUM_TILES_x * $NUM_TILES_y = $TOTAL_TILES tiles to download"
    579592
    580593if [ "$TOTAL_TILES" -gt 200 ] ; then
     
    591604fi
    592605
    593 # end of todo code dump
    594 #######################################################
    595 
    596 
    597 
    598 
    599606
    600607# make a temporary directory
     
    606613rm -f "$TMPDIR"
    607614mkdir "$TMPDIR"
     615WORKDIR=`pwd`
     616cd "$TMPDIR"
    608617
    609618
     
    615624    STYLE="$GIS_OPT_TMBAND"
    616625    TYPE="LandsatTM"
    617     g.message -v "Will download and import $TYPE Data for band $STYLE"
     626    g.message -v "Will download and import $TYPE data for the $STYLE band"
    618627    GetData
    619628    ImportData
     
    624633    STYLE="$GIS_OPT_SRTMBAND"
    625634    TYPE="SRTM"
    626     g.message -v "Will download and import $TYPE Data for band $STYLE"
     635    g.message -v "Will download and import $TYPE data for band $STYLE"
    627636    GetData
    628637    ImportData
     
    633642    STYLE="$GIS_OPT_MONTH"
    634643    TYPE="BMNG"
    635     g.message -v "Will download and import $TYPE Data for the month of $STYLE"
     644    g.message -v "Will download and import $TYPE data for the month of $STYLE"
    636645    GetData
    637646    ImportData
     
    643652    STYLE=""
    644653    TYPE="Daily_Terra"
    645     g.message -v "Will download and import $TYPE Data"
     654    g.message -v "Will download and import $TYPE data"
    646655    GetData
    647656    ImportData
     
    653662    STYLE=""
    654663    TYPE="Daily_Aqua"
    655     g.message -v "Will download and import $TYPE Data"
     664    g.message -v "Will download and import $TYPE data"
    656665    GetData
    657666    ImportData
    658667fi
    659 
     668EXITCODE=$?
    660669
    661670#remove the temp dir
     671cd ..
    662672rm -rf "$TMPDIR"
    663673
     
    665675g.remove region="tmp_tiled_wms.$$" --quiet
    666676
    667 g.message -v "$0 Done."
     677if [ "$EXITCODE" -eq 0 ] ; then
     678   g.message "`basename $0` done."
     679fi
    668680
    669681exit
     682
     683
     684
     685
     686
     687
     688#######################################################
     689# start of old code
     690
     691if [ 0 -eq 1 ] ; then
     692    #There is a bug in nasa WMS service, it provides images which are lager then
     693    #the world :(, we have to crop the images
     694    if [ "$n" = "90" -a "$s" = "-90" ] && \
     695       [ "$w" = "-180" -a "$e" = "180" ] ; then
     696
     697        # check if we have bc
     698        if [ ! -x "`which bc`" ] ; then
     699            g.message -e "bc required, please install first"
     700            exit 1
     701        fi
     702        #We request a smaller image from the wms server
     703        n=`echo "$n - 0.001" | bc`
     704        s=`echo "$s + 0.001" | bc`
     705        e=`echo "$e - 0.001" | bc`
     706        w=`echo "$w + 0.001" | bc`
     707    fi
     708fi
     709
     710# end of old code
     711#######################################################
Note: See TracChangeset for help on using the changeset viewer.