[GRASS-SVN] r52186 - grass-addons/grass6/raster/r.in.onearth
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Jun 22 03:50:45 PDT 2012
Author: hamish
Date: 2012-06-22 03:50:44 -0700 (Fri, 22 Jun 2012)
New Revision: 52186
Added:
grass-addons/grass6/raster/r.in.onearth/r.in.onearth
Removed:
grass-addons/grass6/raster/r.in.onearth/README
grass-addons/grass6/raster/r.in.onearth/r.in.onearth.twms
Modified:
grass-addons/grass6/raster/r.in.onearth/description.html
grass-addons/grass6/raster/r.in.onearth/r.in.onearth.wms
Log:
rename to main script now that it works, cleanup old files a bit
Deleted: grass-addons/grass6/raster/r.in.onearth/README
===================================================================
--- grass-addons/grass6/raster/r.in.onearth/README 2012-06-22 10:47:47 UTC (rev 52185)
+++ grass-addons/grass6/raster/r.in.onearth/README 2012-06-22 10:50:44 UTC (rev 52186)
@@ -1,4 +0,0 @@
-r.in.onearth will need the newest grass6.1-cvs version of g.region. It is in
-CVS since 2005-12-21.
-
-Enjoy
Modified: grass-addons/grass6/raster/r.in.onearth/description.html
===================================================================
--- grass-addons/grass6/raster/r.in.onearth/description.html 2012-06-22 10:47:47 UTC (rev 52185)
+++ grass-addons/grass6/raster/r.in.onearth/description.html 2012-06-22 10:50:44 UTC (rev 52186)
@@ -27,6 +27,7 @@
earth data set including seasonal dynamics.</li>
</ul>
+
<H2>NOTES</H2>
The images will be downloaded and stored into an temporary
@@ -46,6 +47,7 @@
parameters. If these parameters change, the script has to be
modified.
+
<h2>EXAMPLE</h2>
Download Landsat Global Mosaik for Spearfish (SD, USA) area:
@@ -61,6 +63,8 @@
<H2>AUTHOR</H2>
-Soeren Gebbert, Markus Neteler
+Original version by Soeren Gebbert, Markus Neteler<br>
+Rewritten to support tiled WMS by Hamish Bowman
-<p><i>Last changed: $Date$</i>
+<p>
+<i>Last changed: $Date$</i>
Copied: grass-addons/grass6/raster/r.in.onearth/r.in.onearth (from rev 52185, grass-addons/grass6/raster/r.in.onearth/r.in.onearth.twms)
===================================================================
--- grass-addons/grass6/raster/r.in.onearth/r.in.onearth (rev 0)
+++ grass-addons/grass6/raster/r.in.onearth/r.in.onearth 2012-06-22 10:50:44 UTC (rev 52186)
@@ -0,0 +1,711 @@
+#!/bin/sh
+#############################################################################
+# Download and import satellite images direct from the #
+# NASA onearth WMS server into GRASS. #
+# Written by Soeren Gebbert 11/2005 soerengebbert AT gmx de #
+# and Markus Neteler. #
+# Rewritten to support for pre-tiled WMS server by Hamish Bowman #
+# #
+# COPYRIGHT: (C) 2005-2012 by the GRASS Development Team #
+# #
+# This program is free software under the GNU General Public #
+# License (>=v2). Read the file COPYING that comes with GRASS #
+# for details. #
+# #
+#############################################################################
+
+#%Module
+#% description: Download and import satellite images direct from the NASA OnEarth Tiled WMS server into GRASS or to a GeoTIFF image file.
+#%End
+#%option
+#% key: output
+#% gisprompt: new,cell,raster
+#% type: string
+#% description: Output raster map name prefix
+#% required: no
+#%end
+#%option
+#% key: file
+#% gisprompt: new_file,file,file
+#% type: string
+#% description: Output file name prefix
+#% answer: /tmp/test
+#% required: no
+#%end
+#%flag
+#% key: f
+#% description: Do not import to GRASS, create a tiff file instead.
+#%end
+#%flag
+#% key: l
+#% description: Download and Import WMS Global Mosaic, a High Resolution, Global Coverage, Landsat TM 7 mosaic.
+#%end
+#%flag
+#% key: s
+#% description: Download and Import the radar reflectance images produced by the SRTM mission.
+#%end
+#%flag
+#% key: b
+#% description: Download and Import the Blue Marble Next Generation layer, one for each month of the year.
+#%end
+#%flag
+#% key: t
+#% description: Download and Import the composite of data produced by the MODIS Rapid Response System, from data collected yesterday by the MODIS/Terra.
+#%end
+#%flag
+#% key: a
+#% description: Download and Import the composite of data produced by the MODIS Rapid Response System, from data collected yesterday by the MODIS/Aqua.
+#%end
+#%option
+#% key: tmband
+#% type: string
+#% description: NASA Landsat TM bands
+#% options: Red,Green,Blue,IR1,IR2,IR3,ThL,ThH,Pan,visual,pseudo
+#% required: no
+#%end
+#%option
+#% key: srtmband
+#% type: string
+#% description: Radar reflectance bands
+#% options: default,ss1,ss2,ss3,ss4,all
+#% required: no
+#%end
+#%option
+#% key: month
+#% type: string
+#% description: Blue Marble Next Generation layer
+#% options: Jan,Feb,Mar,Apr,May,Jun,Jul,Aug,Sep,Oct,Nov,Dec
+#% required: no
+#%end
+#%option
+#% key: time
+#% type: string
+#% description: The datum of creation for Aqua or Terra satellite images
+#% answer: 2005-3-24
+#% required: no
+#%end
+#%option
+#% key: wgetopt
+#% type: string
+#% description: Options for wget
+#% answer: -c -t 5
+#% required: no
+#%end
+
+
+# Only run if started in GRASS
+if [ -z "$GISBASE" ] ; then
+ echo "You must be in GRASS GIS to run this program." >&2
+ exit 1
+fi
+
+# Parse the arguments
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+
+
+# Set up important vars first
+BASE_URL="http://onearth.jpl.nasa.gov/wms.cgi"
+SRS="EPSG:4326" # This is Lat/Long WGS84
+FORMAT="image/jpeg"
+TILE_SIZE=512
+TYPE="" # Pretty name of layer, no special chars
+TIME="" # &time=YYYY-MM-DD time string for the MODIS time= option
+# Layer IDs:
+GLOBAL_MOSAIC_LAYER="global_mosaic_base" # LANDSAT Global Mosaic
+SRTM_MAG_LAYER="srtm_mag" # srtm_mag (Space Shuttle Radar Topography Mission elevation)
+BMNG_LAYER="BMNG" # NASA' s Blue Marble Next Generation
+DAILY_TERRA_LAYER="daily_terra" # Daily MODIS/TERRA
+DAILY_AQUA_LAYER="daily_aqua" # Daily MODIS/AQUA
+
+IMPORT="true"
+USE_GDALWARP="true" #if we dont have gdalwarp, only LatLong projection is supported
+
+# try a few at once, but don't get abusive
+MAX_DL_JOBS=4
+
+
+# check if we have wget
+if [ ! -x "`which wget`" ] ; then
+ g.message -e "wget required, please install first"
+ exit 1
+fi
+
+# check if we have gdalwarp
+if [ ! -x "`which gdalwarp`" ] ; then
+ g.message -w "gdalwarp is recommended, please install first (script still works in LatLong locations)"
+ USE_GDALWARP="false" # use only LatLong
+fi
+
+
+# Get the data from the NASA server
+GetData()
+{
+ IMPORT="true" #default
+
+ # Compose the filename for the assembled GeoTIFF
+ IMAGEFILE="$TMPDIR/OnEarth_${LAYER}_${STYLE}_$$.tif"
+
+ g.message -v "Download Data ..."
+ g.message -v "Requesting Data from '$BASE_URL'"
+
+ #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"
+ BASE_REQUEST="request=GetMap&layers=$LAYER&srs=$SRS&width=$TILE_SIZE&height=$TILE_SIZE"
+ BASE_META="format=$FORMAT&version=1.1.1&styles=$STYLE$TIME"
+ # what the request would look like if we had a general-purpose WMS server:
+ #BBOX="bbox=$w,$s,$e,$n"
+ #REQ_STRING="$BASE_REQUEST&$BBOX&$BASE_META"
+ #WMS_URL="$BASE_URL?$REQ_STRING"
+
+ #download the File from the Server
+ #wget $WGET_OPTIONS --post-data="$REQ_STRING" "$BASE_URL" -O "$IMAGEFILE"
+
+ for i in `seq $NUM_TILES_x` ; do
+ i_str=`echo "$i" | awk '{printf("%04d", $1)}'`
+
+ # avoid tall+thin jobs from running away
+ MODULUS=`echo "$i $MAX_DL_JOBS" | awk '{print $1 % $2}'`
+ if [ "$MODULUS" = "0" ] ; then
+ wait
+ fi
+
+ for j in `seq $NUM_TILES_y` ; do
+ j_str=`echo "$j" | awk '{printf("%04d", $1)}'`
+ minx=$(echo "$w $i $tile_res" | awk '{printf("%.9f", $1 + (($2 - 1) * $3))}')
+ maxx=$(echo "$w $i $tile_res" | awk '{printf("%.9f", $1 + ($2 * $3))}')
+ maxy=$(echo "$n $j $tile_res" | awk '{printf("%.9f", $1 - (($2 - 1) * $3))}')
+ miny=$(echo "$n $j $tile_res" | awk '{printf("%.9f", $1 - ($2 * $3))}')
+
+ BBOX_TILE="$minx,$miny,$maxx,$maxy"
+ BBOX="bbox=$BBOX_TILE"
+ g.message -d message=" $BBOX"
+
+ REQ_STRING="$BASE_REQUEST&$BBOX&$BASE_META"
+ g.message -d message="URL is [$BASE_URL?$REQ_STRING]"
+
+ #CMD="wget -nv \"$URL\" -O twms_${i_str}_${j_str}.jpg"
+ CMD="wget -nv --post-data='$REQ_STRING' '$BASE_URL' -O 'twms_${i_str}_${j_str}.jpg'"
+
+ MODULUS=`echo "$j $MAX_DL_JOBS" | awk '{print $1 % $2}'`
+ g.message -d message="modulus=$MODULUS"
+ if [ "$MODULUS" = "0" ] ; then
+ # stall to let the background jobs finish
+ g.message "Downloading tile [$i,$j] of [${NUM_TILES_x}x$NUM_TILES_y]."
+ eval $CMD
+ EXITCODE=$?
+ wait
+ if [ "$EXITCODE" -ne 0 ] ; then
+ #retry
+ sleep 1
+ eval $CMD
+ if [ "$EXITCODE" -ne 0 ] ; then
+ g.message -e "wget was not able to download the data"
+ IMPORT="false"
+ return 1
+ fi
+ fi
+ else
+ g.message "Downloading tile [$i,$j] of [${NUM_TILES_x}x$NUM_TILES_y],"
+ eval $CMD &
+ fi
+ done
+ done
+ wait
+
+ #### some checks before going on
+ if [ `ls -1 twms_*.jpg | wc -l` -lt "$TOTAL_TILES" ] ; then
+ g.message -e "Tile(s) appear to be missing."
+ IMPORT="false"
+ return 1
+ fi
+
+ for file in twms_*.jpg ; do
+ if [ ! -s "$file" ] ; then
+ g.message -e "<$file> appears to be empty."
+ IMPORT="false"
+ return 1
+ fi
+ if [ `file "$file" | grep -c JPEG` -eq 0 ] ; then
+ g.message -e "<$file> appears to be bogus."
+ if [ `file "$file" | grep -c XML` -eq 1 ] && [ "$GRASS_VERBOSE" = "true" ] ; then
+ g.message message="--------------------------------------"
+ cat "$file"
+ g.message message="--------------------------------------"
+ fi
+ IMPORT="false"
+ return 1
+ fi
+ done
+
+
+ g.message -v "Converting to pnm ..."
+ for file in twms_*.jpg ; do
+ jpegtopnm "$file" > `basename "$file" .jpg`.pnm
+ done
+
+ g.message -v "Patching ..."
+
+ for j in `seq $NUM_TILES_y` ; do
+ j_str=`echo "$j" | awk '{printf("%04d", $1)}'`
+ tilefiles=`ls twms_*_${j_str}.pnm`
+ #echo $tilefiles
+ pnmcat -lr $tilefiles > "row_${j_str}.pnm"
+ done
+
+ pnmcat -tb row_*.pnm | pnmtojpeg --quality=75 > "mosaic_$$.jpg"
+
+ rm -f twms_*.pnm twms_*.jpg row*.pnm
+
+ #not needed
+ #far_e=$(echo "$e $tile_res" | awk '{printf("%.9f", $1 + $2)}')
+ #far_s=$(echo "$s $tile_res" | awk '{printf("%.9f", $1 - $2)}')
+
+ # Convert to GeoTIFF
+ # gdalwarp can get confused if the file already exists (or at least it used to)
+ if [ -e "mosaic_$$.tif" ] ; then
+ g.message -e "A file called [$TMPDIR/mosaic_$$.tif] already exists. Aborting."
+ exit 1
+ fi
+ gdal_translate "mosaic_$$.jpg" "mosaic_$$.tif" \
+ -a_srs EPSG:4326 -a_ullr $w $n $e $s
+ #-co COMPRESS=DEFLATE
+
+ mv "mosaic_$$.tif" "$IMAGEFILE"
+ \rm "mosaic_$$.jpg"
+
+ if [ -f "$IMAGEFILE" ] ; then
+ IMPORT="true"
+ else
+ g.message -e "wget was not able to download the data"
+ IMPORT="false"
+ return 1
+ fi
+
+ return 0
+}
+
+
+
+#warp the data to the current grass locationa via gdalwarp
+WarpData()
+{
+ if [ "$USE_GDALWARP" = "true" ] ; then
+ g.message -v "Reprojecting the image ..."
+ #create the new imagename
+ IMAGE_WARPED="$TMPDIR/`basename "$IMAGEFILE" .tif`_warped.tif"
+
+ #convert the data to the current location, create Erdas Imagine Images (HFA)
+ gdalwarp -s_srs "$SRS" -t_srs "`g.proj -wf`" \
+ -dstalpha "$IMAGEFILE" "$IMAGE_WARPED"
+ if [ $? -ne 0 ] ; then
+ g.message -i 'Reprojection failed. Aborting.'
+ exitprocedure
+ fi
+ g.message -v "Reprojection successful"
+ #remove the old image and convert the name
+ rm -f "$IMAGEFILE"
+ IMAGEFILE="$IMAGE_WARPED"
+ return 0
+ fi
+ return 1
+}
+
+
+#Import the Data with r.in.gdal
+ImportData()
+{
+ if [ "$IMPORT" != "true" ] ; then
+ return 0
+ fi
+
+ #Check if Tiff file
+ FILETYPE=`file "$IMAGEFILE" | cut -f2 -d':'`
+ g.message -v "File of Type: $FILETYPE"
+ if [ `echo "$FILETYPE" | grep -c 'TIFF'` -ne 1 ] ; then
+ g.message -w "Downloaded file does not appear to be a TIFF image, but will try to import anyway"
+ fi
+
+ g.message -v "Checking Data ..."
+
+ gdalinfo "$IMAGEFILE" > /dev/null
+
+ if [ $? -ne 0 ] ; then
+ g.message -e "Downloaded file is not supported by GDAL, or cannot be imported"
+ fi
+
+ g.message -d "Data checked & was ok."
+
+
+ ### Copy or import
+ if [ "$GIS_FLAG_F" -eq 1 ] ; then
+ # Copy the data to the output file
+ OUTNAME=`echo "$GIS_OPT_OUTPUT" | sed -e 's/_$//'`
+ OUTFILE="${OUTNAME}_${TYPE}_$STYLE.tif"
+
+ g.message -v "Creating output file [$OUTFILE]"
+ cp -f "$IMAGEFILE" "$WORKDIR/$OUTFILE"
+ else
+ # Warp the data!
+ WarpData
+
+ g.message -v "Importing image ..."
+
+ OUTNAME=`basename "$GIS_OPT_OUTPUT" _`
+ OUTNAME="${OUTNAME}_${TYPE}_$STYLE"
+
+ r.in.gdal input="$IMAGEFILE" output="$OUTNAME.tmp_$$" --quiet
+
+ # crop away no data introduced by differing convergence angle
+ g.region "rast=$OUTNAME.tmp_$$.alpha"
+ for COLOR in red green blue ; do
+ r.mapcalc "$OUTNAME.$COLOR = \
+ if($OUTNAME.tmp_$$.alpha == 0, null(), $OUTNAME.tmp_$$.$COLOR)" &
+ done
+ wait
+
+ for COLOR in red green blue ; do
+ r.colors "$OUTNAME.$COLOR" color=grey255 --quiet
+ r.support "$OUTNAME.$COLOR" \
+ title="NASA OnEarth WMS $LAYER ($STYLE)"
+ source1="Data downloaded from NASA's OnEarth WMS server" \
+ source2="$LAYER ($STYLE)" \
+ description="generated by r.in.onearth"
+ done
+
+ g.mremove -f rast="$OUTNAME.tmp_$$.*" --quiet
+ fi
+
+ # free up the disk space ASAP
+ rm -f "$IMAGEFILE"
+ return 0
+}
+
+
+# what to do in case of user break:
+exitprocedure()
+{
+ g.message 'User break!'
+ cd ..
+ rm -rf "$TMPDIR"
+ unset WIND_OVERRIDE
+ g.remove region="tmp_tiled_wms.$$" --quiet
+ exit 1
+}
+trap "exitprocedure" 2 3 15
+
+
+#At least one flag should be set
+if [ $GIS_FLAG_L -eq 0 -a $GIS_FLAG_S -eq 0 -a $GIS_FLAG_B -eq 0 ] \
+ && [ $GIS_FLAG_T -eq 0 -a $GIS_FLAG_A -eq 0 ] ; then
+ g.message -e "Select a flag to specify map type"
+ exit 1
+fi
+
+#Check if a file or a map should be created
+if [ "$GIS_FLAG_F" -eq 1 ] ; then
+ if [ -z "$GIS_OPT_FILE" ] ; then
+ g.message -e "Please specify an output filename"
+ exit 1
+ fi
+ if [ -e "$GIS_OPT_FILE" ] ; then
+ g.message -e "Output filename already exists. Will not overwrite. Aborting."
+ exit 1
+ fi
+fi
+
+#Some mapset informations
+eval `g.gisenv`
+: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
+LOCATION="$GISDBASE/$LOCATION_NAME/$MAPSET"
+PERM="$GISDBASE/$LOCATION_NAME/PERMANENT"
+
+#wget has many options
+WGET_OPTIONS="$GIS_OPT_WGETOPT"
+
+# Tiled WMS is locked at one size fits all
+WIDTH="$TILE_SIZE"
+HEIGHT="$TILE_SIZE"
+
+
+# are we in the server's native projection?
+PROJ_CODE=`g.region -p | grep '^projection:' | awk '{print $2}'`
+DATUM=`g.region -p | grep '^datum:' | awk '{print $2}'`
+if [ $PROJ_CODE -eq 3 ] && [ `echo "$DATUM" | grep -ci 'wgs84'` -eq 1 ] ; then
+ IS_4326="true"
+else
+ IS_4326="false"
+fi
+
+
+#Now get the LatLong Boundingbox
+if [ "$IS_4326" = "true" ] ; then
+ #We have LatLong projection, no warp is needed!
+ USE_GDALWARP="false"
+else
+ eval `g.region -gb`
+ n="$ll_n"
+ s="$ll_s"
+ e="$ll_e"
+ w="$ll_w"
+ g.message -d "Lat-Long WGS84 bounding box was = N $n S $s W $w E $e"
+fi
+
+
+#Break If we have no warp and no LatLong
+if [ "$IS_4326" = "false" ] && [ "$USE_GDALWARP" = "false" ] ; then
+ g.message -e "NASA OnEarth data are in Latitude-Longitude WGS84. The \
+ current location projection differs and you don't \
+ have 'gdalwarp' installed. Aborting."
+ exit 1
+fi
+
+eval `g.region -g | grep 'res='`
+
+if [ "$nsres" != "$ewres" ] ; then
+ g.message -e "East-West and North-South region resolution must be the same"
+ exit 1
+else
+ res="$nsres"
+fi
+
+if [ $PROJ_CODE -ne 3 ] ; then
+ # if not lat/lon, convert map units resolution to degrees
+ TO_METERS=`g.proj -p | grep '^meters.*:' | awk '{print $3}'`
+ res=`echo "$nsres $TO_METERS" | awk '{ printf("%.11f", $1 / ($2 * 1852.0 * 60))}'`
+
+fi
+
+
+#### find appropriate tile resolution
+fixed_res="
+0.000244140625
+0.00048828125
+0.0009765625
+0.001953125
+0.00390625
+0.0078125
+0.015625
+0.03125
+0.0625
+0.125
+0.25
+"
+i=-3
+TAKE=9999
+TAKE_i=9999
+MIN=9999
+for val in $fixed_res ; do
+ # find nearest resolution:
+ #DIFF=`echo "$res $val" | awk '{printf("%.15g", $1 > $2 ? $1 - $2 : $2 - $1)}'`
+ # find nearest finer resolution:
+ DIFF=`echo "$res $val" | awk '{printf("%.15g", $1 > $2 ? $1 - $2 : 9999)}'`
+ NEW_MIN_DIFF=`echo "$MIN $DIFF" | awk '{printf("%d", $1 < $2 ? 0 : 1)}'`
+ if [ "$NEW_MIN_DIFF" -eq 1 ] ; then
+ MIN="$DIFF"
+ TAKE="$val"
+ TAKE_i="$i"
+ fi
+ i=`expr $i + 1`
+done
+if [ "$MIN" = "9999" ] ; then
+ #finer than the finest tested
+ TAKE_i=-3
+ TAKE=0.000244140625
+fi
+tile_res=`echo "$TAKE_i" | awk '{printf("%.15g", 2^$1)}'`
+g.message -d message="min_diff=$MIN [2^$TAKE_i = $tile_res ($TAKE * 512)]"
+
+
+#### adjust region to snap to fixed list of available resolutions
+# setup internal region
+g.region save="tmp_tiled_wms.$$"
+WIND_OVERRIDE="tmp_tiled_wms.$$"
+export WIND_OVERRIDE
+
+
+### now snap the current region to the best-fit tiled grid
+#d.region.box #orig
+if [ $PROJ_CODE -ne 3 ] ; then
+ # if not lat/lon, convert map units resolution to degrees
+ local_res=`echo "$TAKE $TO_METERS" | awk '{ printf("%.15g\n", $1 * $2 * 1852.0 * 60)}'`
+ g.region res="$local_res" -a
+ eval `g.region -gb`
+ n="$ll_n"
+ s="$ll_s"
+ e="$ll_e"
+ w="$ll_w"
+ g.message -d "Lat-Long WGS84 bounding box aligned to = N $n S $s W $w E $e"
+else
+ g.region res="$TAKE" -a
+ eval `g.region -g`
+fi
+#d.redraw
+#d.region.box #after resolution snap
+
+
+
+#### begin snap to nearest tile
+# snap lat,long to upper-left corner of nearest tile grid node
+# Round lat to nearest grid node
+#lat = 90.0 - round( (90.0 - lat) / $res) * $res;
+#lon = -180.0 + floor( (180.0 + lon) / $res) * $res;
+
+# diff is always positive here, so int() acts like floor().
+n=`echo "$n $tile_res" | \
+ awk '{ printf("%.15g", 90.0 - int((90.0 - $1) / $2) * $2) }'`
+s=`echo "$s $tile_res" | awk '
+ function ceil(x)
+ {
+ return (x == int(x)) ? x : int(x)+1
+ }
+ {
+ printf("%.15g", 90.0 - ceil((90.0 - $1) / $2) * $2)
+ }'`
+
+w=`echo "$w $tile_res" | \
+ awk '{ printf("%.15g", -180.0 + int((180.0 + $1) / $2) * $2) }'`
+e=`echo "$e $tile_res" | awk '
+ function ceil(x)
+ {
+ return (x == int(x)) ? x : int(x)+1
+ }
+ {
+ printf("%.15g", -180.0 + ceil((180.0 + $1) / $2) * $2)
+ }'`
+
+g.message -d "Lat-Long WGS84 bounding box snapped to = N $n S $s W $w E $e"
+
+# not needed:
+#g.region n=$n s=$s e=$e w=$w
+#d.redraw
+#d.region.box
+
+
+#NUM_TILES_x=`echo "$e $w $tile_res" | awk '{ printf("%d", 1 + (($1 - $2) / $3)) }'`
+#NUM_TILES_y=`echo "$n $s $tile_res" | awk '{ printf("%d", 1 + (($1 - $2) / $3)) }'`
+NUM_TILES_x=`echo "$e $w $tile_res" | awk '{ printf("%d", ($1 - $2) / $3) }'`
+NUM_TILES_y=`echo "$n $s $tile_res" | awk '{ printf("%d", ($1 - $2) / $3) }'`
+TOTAL_TILES=`expr $NUM_TILES_x \* $NUM_TILES_y`
+
+g.message -v "There are $NUM_TILES_x * $NUM_TILES_y = $TOTAL_TILES tiles to download"
+
+if [ "$TOTAL_TILES" -gt 200 ] ; then
+ g.message -w "You've told it to download $TOTAL_TILES tiles from NASA's server. Are you sure you want to do that?"
+ read result
+ case $result in
+ y | yes | Y | Yes | YES)
+ ;;
+ *)
+ g.message -i "Aborting."
+ exit 1
+ ;;
+ esac
+fi
+
+
+# make a temporary directory
+TMPDIR="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "$TMPDIR" ] ; then
+ g.message -e "Unable to create temporary files"
+ exit 1
+fi
+rm -f "$TMPDIR"
+mkdir "$TMPDIR"
+WORKDIR=`pwd`
+cd "$TMPDIR"
+
+
+# Get the Data and import
+# allow multiple outputs: import every choice that can be made
+
+if [ $GIS_FLAG_L -eq 1 ] ; then
+ LAYER="$GLOBAL_MOSAIC_LAYER"
+ STYLE="$GIS_OPT_TMBAND"
+ TYPE="LandsatTM"
+ g.message -v "Will download and import $TYPE data for the $STYLE band"
+ GetData
+ ImportData
+fi
+
+if [ $GIS_FLAG_S -eq 1 ] ; then
+ LAYER="$SRTM_MAG_LAYER"
+ STYLE="$GIS_OPT_SRTMBAND"
+ TYPE="SRTM"
+ g.message -v "Will download and import $TYPE data for band $STYLE"
+ GetData
+ ImportData
+fi
+
+if [ $GIS_FLAG_B -eq 1 ] ; then
+ LAYER="$BMNG_LAYER"
+ STYLE="$GIS_OPT_MONTH"
+ TYPE="BMNG"
+ g.message -v "Will download and import $TYPE data for the month of $STYLE"
+ GetData
+ ImportData
+fi
+
+if [ $GIS_FLAG_T -eq 1 ] ; then
+ LAYER="$DAILY_TERRA_LAYER"
+ TIME="&time=$GIS_OPT_TIME"
+ STYLE=""
+ TYPE="Daily_Terra"
+ g.message -v "Will download and import $TYPE data"
+ GetData
+ ImportData
+fi
+
+if [ $GIS_FLAG_A -eq 1 ] ; then
+ LAYER="$DAILY_AQUA_LAYER"
+ TIME="&time=$GIS_OPT_TIME"
+ STYLE=""
+ TYPE="Daily_Aqua"
+ g.message -v "Will download and import $TYPE data"
+ GetData
+ ImportData
+fi
+EXITCODE=$?
+
+#remove the temp dir
+cd ..
+rm -rf "$TMPDIR"
+
+unset WIND_OVERRIDE
+g.remove region="tmp_tiled_wms.$$" --quiet
+
+if [ "$EXITCODE" -eq 0 ] ; then
+ g.message "`basename $0` done."
+fi
+
+exit
+
+
+
+
+
+
+#######################################################
+# start of old code
+
+if [ 0 -eq 1 ] ; then
+ #There is a bug in nasa WMS service, it provides images which are lager then
+ #the world :(, we have to crop the images
+ if [ "$n" = "90" -a "$s" = "-90" ] && \
+ [ "$w" = "-180" -a "$e" = "180" ] ; then
+
+ # check if we have bc
+ if [ ! -x "`which bc`" ] ; then
+ g.message -e "bc required, please install first"
+ exit 1
+ fi
+ #We request a smaller image from the wms server
+ n=`echo "$n - 0.001" | bc`
+ s=`echo "$s + 0.001" | bc`
+ e=`echo "$e - 0.001" | bc`
+ w=`echo "$w + 0.001" | bc`
+ fi
+fi
+
+# end of old code
+#######################################################
Deleted: grass-addons/grass6/raster/r.in.onearth/r.in.onearth.twms
===================================================================
--- grass-addons/grass6/raster/r.in.onearth/r.in.onearth.twms 2012-06-22 10:47:47 UTC (rev 52185)
+++ grass-addons/grass6/raster/r.in.onearth/r.in.onearth.twms 2012-06-22 10:50:44 UTC (rev 52186)
@@ -1,711 +0,0 @@
-#!/bin/sh
-#############################################################################
-# Download and import satellite images direct from the #
-# NASA onearth WMS server into GRASS. #
-# Written by Soeren Gebbert 11/2005 soerengebbert AT gmx de #
-# and Markus Neteler. #
-# Rewritten to support for pre-tiled WMS server by Hamish Bowman #
-# #
-# COPYRIGHT: (C) 2005-2012 by the GRASS Development Team #
-# #
-# This program is free software under the GNU General Public #
-# License (>=v2). Read the file COPYING that comes with GRASS #
-# for details. #
-# #
-#############################################################################
-
-#%Module
-#% description: Download and import satellite images direct from the NASA OnEarth Tiled WMS server into GRASS or to a GeoTIFF image file.
-#%End
-#%option
-#% key: output
-#% gisprompt: new,cell,raster
-#% type: string
-#% description: Output raster map name prefix
-#% required: no
-#%end
-#%option
-#% key: file
-#% gisprompt: new_file,file,file
-#% type: string
-#% description: Output file name prefix
-#% answer: /tmp/test
-#% required: no
-#%end
-#%flag
-#% key: f
-#% description: Do not import to GRASS, create a tiff file instead.
-#%end
-#%flag
-#% key: l
-#% description: Download and Import WMS Global Mosaic, a High Resolution, Global Coverage, Landsat TM 7 mosaic.
-#%end
-#%flag
-#% key: s
-#% description: Download and Import the radar reflectance images produced by the SRTM mission.
-#%end
-#%flag
-#% key: b
-#% description: Download and Import the Blue Marble Next Generation layer, one for each month of the year.
-#%end
-#%flag
-#% key: t
-#% description: Download and Import the composite of data produced by the MODIS Rapid Response System, from data collected yesterday by the MODIS/Terra.
-#%end
-#%flag
-#% key: a
-#% description: Download and Import the composite of data produced by the MODIS Rapid Response System, from data collected yesterday by the MODIS/Aqua.
-#%end
-#%option
-#% key: tmband
-#% type: string
-#% description: NASA Landsat TM bands
-#% options: Red,Green,Blue,IR1,IR2,IR3,ThL,ThH,Pan,visual,pseudo
-#% required: no
-#%end
-#%option
-#% key: srtmband
-#% type: string
-#% description: Radar reflectance bands
-#% options: default,ss1,ss2,ss3,ss4,all
-#% required: no
-#%end
-#%option
-#% key: month
-#% type: string
-#% description: Blue Marble Next Generation layer
-#% options: Jan,Feb,Mar,Apr,May,Jun,Jul,Aug,Sep,Oct,Nov,Dec
-#% required: no
-#%end
-#%option
-#% key: time
-#% type: string
-#% description: The datum of creation for Aqua or Terra satellite images
-#% answer: 2005-3-24
-#% required: no
-#%end
-#%option
-#% key: wgetopt
-#% type: string
-#% description: Options for wget
-#% answer: -c -t 5
-#% required: no
-#%end
-
-
-# Only run if started in GRASS
-if [ -z "$GISBASE" ] ; then
- echo "You must be in GRASS GIS to run this program." >&2
- exit 1
-fi
-
-# Parse the arguments
-if [ "$1" != "@ARGS_PARSED@" ] ; then
- exec g.parser "$0" "$@"
-fi
-
-
-# Set up important vars first
-BASE_URL="http://onearth.jpl.nasa.gov/wms.cgi"
-SRS="EPSG:4326" # This is Lat/Long WGS84
-FORMAT="image/jpeg"
-TILE_SIZE=512
-TYPE="" # Pretty name of layer, no special chars
-TIME="" # &time=YYYY-MM-DD time string for the MODIS time= option
-# Layer IDs:
-GLOBAL_MOSAIC_LAYER="global_mosaic_base" # LANDSAT Global Mosaic
-SRTM_MAG_LAYER="srtm_mag" # srtm_mag (Space Shuttle Radar Topography Mission elevation)
-BMNG_LAYER="BMNG" # NASA' s Blue Marble Next Generation
-DAILY_TERRA_LAYER="daily_terra" # Daily MODIS/TERRA
-DAILY_AQUA_LAYER="daily_aqua" # Daily MODIS/AQUA
-
-IMPORT="true"
-USE_GDALWARP="true" #if we dont have gdalwarp, only LatLong projection is supported
-
-# try a few at once, but don't get abusive
-MAX_DL_JOBS=4
-
-
-# check if we have wget
-if [ ! -x "`which wget`" ] ; then
- g.message -e "wget required, please install first"
- exit 1
-fi
-
-# check if we have gdalwarp
-if [ ! -x "`which gdalwarp`" ] ; then
- g.message -w "gdalwarp is recommended, please install first (script still works in LatLong locations)"
- USE_GDALWARP="false" # use only LatLong
-fi
-
-
-# Get the data from the NASA server
-GetData()
-{
- IMPORT="true" #default
-
- # Compose the filename for the assembled GeoTIFF
- IMAGEFILE="$TMPDIR/OnEarth_${LAYER}_${STYLE}_$$.tif"
-
- g.message -v "Download Data ..."
- g.message -v "Requesting Data from '$BASE_URL'"
-
- #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"
- BASE_REQUEST="request=GetMap&layers=$LAYER&srs=$SRS&width=$TILE_SIZE&height=$TILE_SIZE"
- BASE_META="format=$FORMAT&version=1.1.1&styles=$STYLE$TIME"
- # what the request would look like if we had a general-purpose WMS server:
- #BBOX="bbox=$w,$s,$e,$n"
- #REQ_STRING="$BASE_REQUEST&$BBOX&$BASE_META"
- #WMS_URL="$BASE_URL?$REQ_STRING"
-
- #download the File from the Server
- #wget $WGET_OPTIONS --post-data="$REQ_STRING" "$BASE_URL" -O "$IMAGEFILE"
-
- for i in `seq $NUM_TILES_x` ; do
- i_str=`echo "$i" | awk '{printf("%04d", $1)}'`
-
- # avoid tall+thin jobs from running away
- MODULUS=`echo "$i $MAX_DL_JOBS" | awk '{print $1 % $2}'`
- if [ "$MODULUS" = "0" ] ; then
- wait
- fi
-
- for j in `seq $NUM_TILES_y` ; do
- j_str=`echo "$j" | awk '{printf("%04d", $1)}'`
- minx=$(echo "$w $i $tile_res" | awk '{printf("%.9f", $1 + (($2 - 1) * $3))}')
- maxx=$(echo "$w $i $tile_res" | awk '{printf("%.9f", $1 + ($2 * $3))}')
- maxy=$(echo "$n $j $tile_res" | awk '{printf("%.9f", $1 - (($2 - 1) * $3))}')
- miny=$(echo "$n $j $tile_res" | awk '{printf("%.9f", $1 - ($2 * $3))}')
-
- BBOX_TILE="$minx,$miny,$maxx,$maxy"
- BBOX="bbox=$BBOX_TILE"
- g.message -d message=" $BBOX"
-
- REQ_STRING="$BASE_REQUEST&$BBOX&$BASE_META"
- g.message -d message="URL is [$BASE_URL?$REQ_STRING]"
-
- #CMD="wget -nv \"$URL\" -O twms_${i_str}_${j_str}.jpg"
- CMD="wget -nv --post-data='$REQ_STRING' '$BASE_URL' -O 'twms_${i_str}_${j_str}.jpg'"
-
- MODULUS=`echo "$j $MAX_DL_JOBS" | awk '{print $1 % $2}'`
- g.message -d message="modulus=$MODULUS"
- if [ "$MODULUS" = "0" ] ; then
- # stall to let the background jobs finish
- g.message "Downloading tile [$i,$j] of [${NUM_TILES_x}x$NUM_TILES_y]."
- eval $CMD
- EXITCODE=$?
- wait
- if [ "$EXITCODE" -ne 0 ] ; then
- #retry
- sleep 1
- eval $CMD
- if [ "$EXITCODE" -ne 0 ] ; then
- g.message -e "wget was not able to download the data"
- IMPORT="false"
- return 1
- fi
- fi
- else
- g.message "Downloading tile [$i,$j] of [${NUM_TILES_x}x$NUM_TILES_y],"
- eval $CMD &
- fi
- done
- done
- wait
-
- #### some checks before going on
- if [ `ls -1 twms_*.jpg | wc -l` -lt "$TOTAL_TILES" ] ; then
- g.message -e "Tile(s) appear to be missing."
- IMPORT="false"
- return 1
- fi
-
- for file in twms_*.jpg ; do
- if [ ! -s "$file" ] ; then
- g.message -e "<$file> appears to be empty."
- IMPORT="false"
- return 1
- fi
- if [ `file "$file" | grep -c JPEG` -eq 0 ] ; then
- g.message -e "<$file> appears to be bogus."
- if [ `file "$file" | grep -c XML` -eq 1 ] && [ "$GRASS_VERBOSE" = "true" ] ; then
- g.message message="--------------------------------------"
- cat "$file"
- g.message message="--------------------------------------"
- fi
- IMPORT="false"
- return 1
- fi
- done
-
-
- g.message -v "Converting to pnm ..."
- for file in twms_*.jpg ; do
- jpegtopnm "$file" > `basename "$file" .jpg`.pnm
- done
-
- g.message -v "Patching ..."
-
- for j in `seq $NUM_TILES_y` ; do
- j_str=`echo "$j" | awk '{printf("%04d", $1)}'`
- tilefiles=`ls twms_*_${j_str}.pnm`
- #echo $tilefiles
- pnmcat -lr $tilefiles > "row_${j_str}.pnm"
- done
-
- pnmcat -tb row_*.pnm | pnmtojpeg --quality=75 > "mosaic_$$.jpg"
-
- rm -f twms_*.pnm twms_*.jpg row*.pnm
-
- #not needed
- #far_e=$(echo "$e $tile_res" | awk '{printf("%.9f", $1 + $2)}')
- #far_s=$(echo "$s $tile_res" | awk '{printf("%.9f", $1 - $2)}')
-
- # Convert to GeoTIFF
- # gdalwarp can get confused if the file already exists (or at least it used to)
- if [ -e "mosaic_$$.tif" ] ; then
- g.message -e "A file called [$TMPDIR/mosaic_$$.tif] already exists. Aborting."
- exit 1
- fi
- gdal_translate "mosaic_$$.jpg" "mosaic_$$.tif" \
- -a_srs EPSG:4326 -a_ullr $w $n $e $s
- #-co COMPRESS=DEFLATE
-
- mv "mosaic_$$.tif" "$IMAGEFILE"
- \rm "mosaic_$$.jpg"
-
- if [ -f "$IMAGEFILE" ] ; then
- IMPORT="true"
- else
- g.message -e "wget was not able to download the data"
- IMPORT="false"
- return 1
- fi
-
- return 0
-}
-
-
-
-#warp the data to the current grass locationa via gdalwarp
-WarpData()
-{
- if [ "$USE_GDALWARP" = "true" ] ; then
- g.message -v "Reprojecting the image ..."
- #create the new imagename
- IMAGE_WARPED="$TMPDIR/`basename "$IMAGEFILE" .tif`_warped.tif"
-
- #convert the data to the current location, create Erdas Imagine Images (HFA)
- gdalwarp -s_srs "$SRS" -t_srs "`g.proj -wf`" \
- -dstalpha "$IMAGEFILE" "$IMAGE_WARPED"
- if [ $? -ne 0 ] ; then
- g.message -i 'Reprojection failed. Aborting.'
- exitprocedure
- fi
- g.message -v "Reprojection successful"
- #remove the old image and convert the name
- rm -f "$IMAGEFILE"
- IMAGEFILE="$IMAGE_WARPED"
- return 0
- fi
- return 1
-}
-
-
-#Import the Data with r.in.gdal
-ImportData()
-{
- if [ "$IMPORT" != "true" ] ; then
- return 0
- fi
-
- #Check if Tiff file
- FILETYPE=`file "$IMAGEFILE" | cut -f2 -d':'`
- g.message -v "File of Type: $FILETYPE"
- if [ `echo "$FILETYPE" | grep -c 'TIFF'` -ne 1 ] ; then
- g.message -w "Downloaded file does not appear to be a TIFF image, but will try to import anyway"
- fi
-
- g.message -v "Checking Data ..."
-
- gdalinfo "$IMAGEFILE" > /dev/null
-
- if [ $? -ne 0 ] ; then
- g.message -e "Downloaded file is not supported by GDAL, or cannot be imported"
- fi
-
- g.message -d "Data checked & was ok."
-
-
- ### Copy or import
- if [ "$GIS_FLAG_F" -eq 1 ] ; then
- # Copy the data to the output file
- OUTNAME=`echo "$GIS_OPT_OUTPUT" | sed -e 's/_$//'`
- OUTFILE="${OUTNAME}_${TYPE}_$STYLE.tif"
-
- g.message -v "Creating output file [$OUTFILE]"
- cp -f "$IMAGEFILE" "$WORKDIR/$OUTFILE"
- else
- # Warp the data!
- WarpData
-
- g.message -v "Importing image ..."
-
- OUTNAME=`basename "$GIS_OPT_OUTPUT" _`
- OUTNAME="${OUTNAME}_${TYPE}_$STYLE"
-
- r.in.gdal input="$IMAGEFILE" output="$OUTNAME.tmp_$$" --quiet
-
- # crop away no data introduced by differing convergence angle
- g.region "rast=$OUTNAME.tmp_$$.alpha"
- for COLOR in red green blue ; do
- r.mapcalc "$OUTNAME.$COLOR = \
- if($OUTNAME.tmp_$$.alpha == 0, null(), $OUTNAME.tmp_$$.$COLOR)" &
- done
- wait
-
- for COLOR in red green blue ; do
- r.colors "$OUTNAME.$COLOR" color=grey255 --quiet
- r.support "$OUTNAME.$COLOR" \
- title="NASA OnEarth WMS $LAYER ($STYLE)"
- source1="Data downloaded from NASA's OnEarth WMS server" \
- source2="$LAYER ($STYLE)" \
- description="generated by r.in.onearth"
- done
-
- g.mremove -f rast="$OUTNAME.tmp_$$.*" --quiet
- fi
-
- # free up the disk space ASAP
- rm -f "$IMAGEFILE"
- return 0
-}
-
-
-# what to do in case of user break:
-exitprocedure()
-{
- g.message 'User break!'
- cd ..
- rm -rf "$TMPDIR"
- unset WIND_OVERRIDE
- g.remove region="tmp_tiled_wms.$$" --quiet
- exit 1
-}
-trap "exitprocedure" 2 3 15
-
-
-#At least one flag should be set
-if [ $GIS_FLAG_L -eq 0 -a $GIS_FLAG_S -eq 0 -a $GIS_FLAG_B -eq 0 ] \
- && [ $GIS_FLAG_T -eq 0 -a $GIS_FLAG_A -eq 0 ] ; then
- g.message -e "Select a flag to specify map type"
- exit 1
-fi
-
-#Check if a file or a map should be created
-if [ "$GIS_FLAG_F" -eq 1 ] ; then
- if [ -z "$GIS_OPT_FILE" ] ; then
- g.message -e "Please specify an output filename"
- exit 1
- fi
- if [ -e "$GIS_OPT_FILE" ] ; then
- g.message -e "Output filename already exists. Will not overwrite. Aborting."
- exit 1
- fi
-fi
-
-#Some mapset informations
-eval `g.gisenv`
-: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
-LOCATION="$GISDBASE/$LOCATION_NAME/$MAPSET"
-PERM="$GISDBASE/$LOCATION_NAME/PERMANENT"
-
-#wget has many options
-WGET_OPTIONS="$GIS_OPT_WGETOPT"
-
-# Tiled WMS is locked at one size fits all
-WIDTH="$TILE_SIZE"
-HEIGHT="$TILE_SIZE"
-
-
-# are we in the server's native projection?
-PROJ_CODE=`g.region -p | grep '^projection:' | awk '{print $2}'`
-DATUM=`g.region -p | grep '^datum:' | awk '{print $2}'`
-if [ $PROJ_CODE -eq 3 ] && [ `echo "$DATUM" | grep -ci 'wgs84'` -eq 1 ] ; then
- IS_4326="true"
-else
- IS_4326="false"
-fi
-
-
-#Now get the LatLong Boundingbox
-if [ "$IS_4326" = "true" ] ; then
- #We have LatLong projection, no warp is needed!
- USE_GDALWARP="false"
-else
- eval `g.region -gb`
- n="$ll_n"
- s="$ll_s"
- e="$ll_e"
- w="$ll_w"
- g.message -d "Lat-Long WGS84 bounding box was = N $n S $s W $w E $e"
-fi
-
-
-#Break If we have no warp and no LatLong
-if [ "$IS_4326" = "false" ] && [ "$USE_GDALWARP" = "false" ] ; then
- g.message -e "NASA OnEarth data are in Latitude-Longitude WGS84. The \
- current location projection differs and you don't \
- have 'gdalwarp' installed. Aborting."
- exit 1
-fi
-
-eval `g.region -g | grep 'res='`
-
-if [ "$nsres" != "$ewres" ] ; then
- g.message -e "East-West and North-South region resolution must be the same"
- exit 1
-else
- res="$nsres"
-fi
-
-if [ $PROJ_CODE -ne 3 ] ; then
- # if not lat/lon, convert map units resolution to degrees
- TO_METERS=`g.proj -p | grep '^meters.*:' | awk '{print $3}'`
- res=`echo "$nsres $TO_METERS" | awk '{ printf("%.11f", $1 / ($2 * 1852.0 * 60))}'`
-
-fi
-
-
-#### find appropriate tile resolution
-fixed_res="
-0.000244140625
-0.00048828125
-0.0009765625
-0.001953125
-0.00390625
-0.0078125
-0.015625
-0.03125
-0.0625
-0.125
-0.25
-"
-i=-3
-TAKE=9999
-TAKE_i=9999
-MIN=9999
-for val in $fixed_res ; do
- # find nearest resolution:
- #DIFF=`echo "$res $val" | awk '{printf("%.15g", $1 > $2 ? $1 - $2 : $2 - $1)}'`
- # find nearest finer resolution:
- DIFF=`echo "$res $val" | awk '{printf("%.15g", $1 > $2 ? $1 - $2 : 9999)}'`
- NEW_MIN_DIFF=`echo "$MIN $DIFF" | awk '{printf("%d", $1 < $2 ? 0 : 1)}'`
- if [ "$NEW_MIN_DIFF" -eq 1 ] ; then
- MIN="$DIFF"
- TAKE="$val"
- TAKE_i="$i"
- fi
- i=`expr $i + 1`
-done
-if [ "$MIN" = "9999" ] ; then
- #finer than the finest tested
- TAKE_i=-3
- TAKE=0.000244140625
-fi
-tile_res=`echo "$TAKE_i" | awk '{printf("%.15g", 2^$1)}'`
-g.message -d message="min_diff=$MIN [2^$TAKE_i = $tile_res ($TAKE * 512)]"
-
-
-#### adjust region to snap to fixed list of available resolutions
-# setup internal region
-g.region save="tmp_tiled_wms.$$"
-WIND_OVERRIDE="tmp_tiled_wms.$$"
-export WIND_OVERRIDE
-
-
-### now snap the current region to the best-fit tiled grid
-#d.region.box #orig
-if [ $PROJ_CODE -ne 3 ] ; then
- # if not lat/lon, convert map units resolution to degrees
- local_res=`echo "$TAKE $TO_METERS" | awk '{ printf("%.15g\n", $1 * $2 * 1852.0 * 60)}'`
- g.region res="$local_res" -a
- eval `g.region -gb`
- n="$ll_n"
- s="$ll_s"
- e="$ll_e"
- w="$ll_w"
- g.message -d "Lat-Long WGS84 bounding box aligned to = N $n S $s W $w E $e"
-else
- g.region res="$TAKE" -a
- eval `g.region -g`
-fi
-#d.redraw
-#d.region.box #after resolution snap
-
-
-
-#### begin snap to nearest tile
-# snap lat,long to upper-left corner of nearest tile grid node
-# Round lat to nearest grid node
-#lat = 90.0 - round( (90.0 - lat) / $res) * $res;
-#lon = -180.0 + floor( (180.0 + lon) / $res) * $res;
-
-# diff is always positive here, so int() acts like floor().
-n=`echo "$n $tile_res" | \
- awk '{ printf("%.15g", 90.0 - int((90.0 - $1) / $2) * $2) }'`
-s=`echo "$s $tile_res" | awk '
- function ceil(x)
- {
- return (x == int(x)) ? x : int(x)+1
- }
- {
- printf("%.15g", 90.0 - ceil((90.0 - $1) / $2) * $2)
- }'`
-
-w=`echo "$w $tile_res" | \
- awk '{ printf("%.15g", -180.0 + int((180.0 + $1) / $2) * $2) }'`
-e=`echo "$e $tile_res" | awk '
- function ceil(x)
- {
- return (x == int(x)) ? x : int(x)+1
- }
- {
- printf("%.15g", -180.0 + ceil((180.0 + $1) / $2) * $2)
- }'`
-
-g.message -d "Lat-Long WGS84 bounding box snapped to = N $n S $s W $w E $e"
-
-# not needed:
-#g.region n=$n s=$s e=$e w=$w
-#d.redraw
-#d.region.box
-
-
-#NUM_TILES_x=`echo "$e $w $tile_res" | awk '{ printf("%d", 1 + (($1 - $2) / $3)) }'`
-#NUM_TILES_y=`echo "$n $s $tile_res" | awk '{ printf("%d", 1 + (($1 - $2) / $3)) }'`
-NUM_TILES_x=`echo "$e $w $tile_res" | awk '{ printf("%d", ($1 - $2) / $3) }'`
-NUM_TILES_y=`echo "$n $s $tile_res" | awk '{ printf("%d", ($1 - $2) / $3) }'`
-TOTAL_TILES=`expr $NUM_TILES_x \* $NUM_TILES_y`
-
-g.message -v "There are $NUM_TILES_x * $NUM_TILES_y = $TOTAL_TILES tiles to download"
-
-if [ "$TOTAL_TILES" -gt 200 ] ; then
- g.message -w "You've told it to download $TOTAL_TILES tiles from NASA's server. Are you sure you want to do that?"
- read result
- case $result in
- y | yes | Y | Yes | YES)
- ;;
- *)
- g.message -i "Aborting."
- exit 1
- ;;
- esac
-fi
-
-
-# make a temporary directory
-TMPDIR="`g.tempfile pid=$$`"
-if [ $? -ne 0 ] || [ -z "$TMPDIR" ] ; then
- g.message -e "Unable to create temporary files"
- exit 1
-fi
-rm -f "$TMPDIR"
-mkdir "$TMPDIR"
-WORKDIR=`pwd`
-cd "$TMPDIR"
-
-
-# Get the Data and import
-# allow multiple outputs: import every choice that can be made
-
-if [ $GIS_FLAG_L -eq 1 ] ; then
- LAYER="$GLOBAL_MOSAIC_LAYER"
- STYLE="$GIS_OPT_TMBAND"
- TYPE="LandsatTM"
- g.message -v "Will download and import $TYPE data for the $STYLE band"
- GetData
- ImportData
-fi
-
-if [ $GIS_FLAG_S -eq 1 ] ; then
- LAYER="$SRTM_MAG_LAYER"
- STYLE="$GIS_OPT_SRTMBAND"
- TYPE="SRTM"
- g.message -v "Will download and import $TYPE data for band $STYLE"
- GetData
- ImportData
-fi
-
-if [ $GIS_FLAG_B -eq 1 ] ; then
- LAYER="$BMNG_LAYER"
- STYLE="$GIS_OPT_MONTH"
- TYPE="BMNG"
- g.message -v "Will download and import $TYPE data for the month of $STYLE"
- GetData
- ImportData
-fi
-
-if [ $GIS_FLAG_T -eq 1 ] ; then
- LAYER="$DAILY_TERRA_LAYER"
- TIME="&time=$GIS_OPT_TIME"
- STYLE=""
- TYPE="Daily_Terra"
- g.message -v "Will download and import $TYPE data"
- GetData
- ImportData
-fi
-
-if [ $GIS_FLAG_A -eq 1 ] ; then
- LAYER="$DAILY_AQUA_LAYER"
- TIME="&time=$GIS_OPT_TIME"
- STYLE=""
- TYPE="Daily_Aqua"
- g.message -v "Will download and import $TYPE data"
- GetData
- ImportData
-fi
-EXITCODE=$?
-
-#remove the temp dir
-cd ..
-rm -rf "$TMPDIR"
-
-unset WIND_OVERRIDE
-g.remove region="tmp_tiled_wms.$$" --quiet
-
-if [ "$EXITCODE" -eq 0 ] ; then
- g.message "`basename $0` done."
-fi
-
-exit
-
-
-
-
-
-
-#######################################################
-# start of old code
-
-if [ 0 -eq 1 ] ; then
- #There is a bug in nasa WMS service, it provides images which are lager then
- #the world :(, we have to crop the images
- if [ "$n" = "90" -a "$s" = "-90" ] && \
- [ "$w" = "-180" -a "$e" = "180" ] ; then
-
- # check if we have bc
- if [ ! -x "`which bc`" ] ; then
- g.message -e "bc required, please install first"
- exit 1
- fi
- #We request a smaller image from the wms server
- n=`echo "$n - 0.001" | bc`
- s=`echo "$s + 0.001" | bc`
- e=`echo "$e - 0.001" | bc`
- w=`echo "$w + 0.001" | bc`
- fi
-fi
-
-# end of old code
-#######################################################
Modified: grass-addons/grass6/raster/r.in.onearth/r.in.onearth.wms
===================================================================
--- grass-addons/grass6/raster/r.in.onearth/r.in.onearth.wms 2012-06-22 10:47:47 UTC (rev 52185)
+++ grass-addons/grass6/raster/r.in.onearth/r.in.onearth.wms 2012-06-22 10:50:44 UTC (rev 52186)
@@ -1,4 +1,9 @@
#!/bin/sh
+#
+# The OnEarth server changed, it no longer supports arbitrary WMS requests.
+# So this script doesn't work any more, but is a nice example of how to
+# interface with a WMS server so remains here for educational purposes.
+#
#############################################################################
# Download and import satellite images direct from the #
# NASA onearth WMS server into GRASS. #
@@ -12,6 +17,10 @@
# for details. #
# #
#############################################################################
+# r.in.onearth will need the newest grass6.1-cvs version of g.region. It is in
+# CVS since 2005-12-21.
+#
+# Enjoy
#%Module
#% description: Download and import satellite images direct from the NASA onearth WMS server into GRASS or to a geo-tiff image file.
More information about the grass-commit
mailing list