Re: adding illumination to hemisphere (-JG)

gmt_shell_functions.sh is part of the GMT distribution.

On Mar 26, 2014, at 6:51 AM, Andy Jacobson <[email protected]> wrote:

> Hi Keith,
> 
> Thanks for showing us this script.  I would like to try it out, but it appears that we may also need your gmt_shell_functions.sh file.  Are there interesting nuggets of wisdom in that file, too?
> 
> Best,
> 
> Andy
> 
> On Wed 26 Mar 2014, at 09:16 , keithpickering . <[email protected]> wrote:
> 
>> Hi Pasquale,
>> 
>> I created a script to do that some time ago, which I include below. The script depends on the "bc" utility to handle the math, but most Unix/Linux installations should have that.
>> 
>> The script actually loops through various viewpoints and subsolar points to create frames that form a "spinning earth" .gif, but you can easily modify it to choose just one frame or view.
>> 
>> Best wishes, 
>> 
>> Keith Pickering
>> 
>> *****************************
>> 
>> #!/bin/sh
>> ##############################################################
>> #            Animated turning globe with terminator
>> #            Author: Keith Pickering
>> #
>> ##############################################################
>> # Usage: set latitude and longitude of subsolar point, plus
>> #        latitude and longitude of viewpoint.
>> # Longitudes of the two must be marginally different, or an
>> # error is thrown. Also, direct 180 opposite longitudes 
>> # throws the same error. Viewpoints over North and South poles
>> # give warnings, but produce good output. All other rotations
>> # should work with no errors or warnings.
>> ###############################################################
>> 
>> # Create file name in tif format
>> height=1.0i
>> width=1.0i
>> dpi=100
>> . gmt_shell_functions.sh
>> gmtset DOTS_PR_INCH $dpi PAPER_MEDIA=Custom_${width}x${height}
>> name=globe3
>> # Set starting frame number, total number of frames, and number of degrees of rotation between each frame
>> frame=1
>> nframes=90
>> framesd=4  
>> mkdir $$
>> # Standard variables, non-looping:
>> # Define latitide and longitude of subsolar point:
>> slon=164.2152
>> slat=15.6865
>> # Define latitude and longitude of view center:
>> vlon=-83
>> vlat=41
>> # Mapping options and colors
>> D="-Dc -A600/0/1"
>> day="-G110/150/120 -S33/33/150"
>> nite="-G33/45/36 -S9/9/45"
>> dr=.017453293
>> #
>> # Main loop starts here:
>> #
>> while [ $frame -le $nframes ]; do
>>   ps=$$/$(eval "gmt_set_framename $name $frame").ps
>>   vlon=$(echo "-81.5-($frame*$framesd) "|bc)
>>   slon=$(echo "sl=$vlon+241.2; if(sl<-180) sl+=360;sl"|bc)
>>   echo "file: " $ps "vlon: " $vlon "slon: " $slon
>> OPT="-Rd -JG"$vlon/$vlat"/.9i -P"
>> 
>> ## Standard terminator script starts here
>> tlat=$(echo "scale=4; $slat-90"|bc) 
>> vlaty=$(echo "scale=4; $vlat-90"|bc)
>> ld=$(bc << EOF
>> scale=4
>> ld=$slon-($vlon)
>> if (ld>180) ld-=360
>> if (ld<-180) ld+=360
>> ld
>> EOF
>> )
>> # Gibbous/crescent switch: 1=crescent
>> # New routine: compute spherical angle between viewpoint and subsolar point.
>> ldi=$(bc -l << EOF
>> scale=4
>> ld=$slon-($vlon)
>> cd=s($dr*$vlat)*s($dr*$slat)+c($dr*$vlat)*c($dr*$slat)*c($dr*ld) 
>> if (cd < 0) ldi=1 else ldi=0
>> ldi
>> EOF
>> )
>> # Angle for -L option in clip path.
>> l1=$(bc -l << EOF
>> scale=4
>> sb=c($slat*$dr)*c($vlat*$dr)*c($ld*$dr)+s($slat*$dr)*s($vlat*$dr)
>> cb=sqrt(1-sb*sb)
>> b=a(sb/cb)
>> y=s($slat*$dr)-sb*s($vlat*$dr)
>> x=c($slat*$dr)*s($ld*$dr)*c($vlat*$dr)
>> res=a(y/x)
>> # resolve arctan ambiguity:
>> if (x<0) if (y<0) res-=(180*$dr) else res+=(180*$dr)
>> res/$dr
>> EOF
>> )
>> 
>> # generate terminator clipping path
>> project -G5 -S -T$slon/$slat -C$slon/$tlat -L0/360 > clip1.txt
>> 
>> if [ $ldi -eq 1 ]; then
>> # cresent: central meridian is unlit
>>   l2=$(echo "scale=4; $l1-180"|bc -l) 
>>   # GC2. At right angles to the view, covering the lit half of the globe.
>>   project -G5 -S -T$vlon/$vlat -C$vlon/$vlaty -L$l2/$l1 > clip2.txt
>>   # day hemisphere: whole globe, will be mostly overlain with night color scheme.
>>   pscoast $OPT $day $D -K > $ps
>>   # Terminator ellipse
>>   psclip clip1.txt  $OPT -K -O >> $ps
>>   pscoast  $OPT $nite $D -K -O >> $ps
>>   psclip -C -O -K >> $ps 
>>   # Unlit half.
>>   psclip clip2.txt $OPT -K -O >> $ps
>>   pscoast  $OPT $nite $D -K -O >> $ps
>>   psclip -C -O >> $ps
>> else
>> # gibbous: central meridian is lit
>>   l2=$(echo "scale=4; $l1+180"|bc -l) 
>>   # GC2. At right angles to the view, covering the lit half of the globe.
>>   project -G5 -S -T$vlon/$vlat -C$vlon/$vlaty -L$l1/$l2 > clip2.txt
>>   # night hemisphere: whole globe, will be mostly overlain with day color scheme.
>>   pscoast $OPT $nite $D -K > $ps
>>   # Terminator ellipse
>>   psclip clip1.txt  $OPT -K -O >> $ps
>>   pscoast  $OPT $day $D -K -O >> $ps
>>   psclip -C -O -K >> $ps 
>>   # Lit half.
>>   psclip clip2.txt $OPT -K -O >> $ps
>>   pscoast  $OPT $day $D -K -O >> $ps
>>   psclip -C -O >>$ps
>> fi
>>   ps2raster -E$dpi -Tt $ps
>>   frame=$(eval "gmt_set_framenext $frame")
>> done
>> echo "converting to gif, this may take a while ..."
>>   convert -delay 4 -loop 0 $$/*.tif $name.gif
>> echo $name".gif conversion, done, cleaning up ..."
>>   rm $$/*.ps
>> cat << EOF > $name.html
>> <HTML>
>> <TITLE>Rotating Globe</TITLE>
>> <BODY>
>> <img src="$name.gif" border="1">
>> </BODY>
>> </HTML>
>> EOF
>> 
>> 
>> On Tue, Mar 25, 2014 at 8:58 PM, Pasquale Tricarico <[email protected]> wrote:
>> Hi,
>> 
>> I am using GMT 4.5.12, using grdimage and -JG projection to obtain an
>> hemisphere view of a planet. I would like to add illumination, i.e. by
>> specifying
>> the lat/lon of the sun, so that the correct part of the hemisphere is
>> illuminated and the other part is much darker, but still visible. It
>> looks like grdgradient may not be able to do this as it would have a
>> fixed azimuth for illumination, maybe by using grdmath? Any
>> suggestions will be greatly appreciated.
>> 
>> Regards,
>> Pasquale
>> 
>> --
>> Pasquale Tricarico, Ph.D.
>> Planetary Science Institute
>> http://orbit.psi.edu/~tricaric
>> 
>> To unsubscribe, send the message "signoff gmt-help" to [email protected]
>> 
>> To unsubscribe, send the message "signoff gmt-help" to [email protected]
> 
> -- 
> Andy Jacobson
> [email protected]
> 
> NOAA Earth System Research Lab
> Global Monitoring Division
> 325 Broadway R/GMD1
> Boulder, Colorado 80305
> 
> 303/578-2237
> 
> To unsubscribe, send the message "signoff gmt-help" to [email protected]

To unsubscribe, send the message "signoff gmt-help" to [email protected]