Earth at night

From GRASS-Wiki
Revision as of 12:50, 15 September 2011 by ⚠️AnneGhisla (talk | contribs) (→‎Get data: fixed broken link)
Jump to navigation Jump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.

Import NASA's Night on Earth and Reproject in Polar Stereographic

This micro-tutorial requires GRASS 6.5

Get data

Download http://eoimages.gsfc.nasa.gov/ve/1438/land_ocean_ice_lights_2048.tif from NASA's Visible Earth website.

Import

  • Create a WGS84 Lat/Lon location (EPSG # 4326)
# import into lat/lon, overriding projection check and
# setting dummy region bounds which fit into geographic
# space (the file is missing georeferencing metadata but
# we do know exactly what it should be).

r.in.gdal -o -l in=land_ocean_ice_lights_8192.tif out=land_ocean_ice_lights

for BAND in red green blue ; do
  r.region map=land_ocean_ice_lights.$BAND n=90 s=-90 w=-180 e=180
done

g.region n=90 s=45 w=-180 e=180

# create 15° and 30° grid line vectors
v.mkgrid map=grid_15 box=15,15 grid=12,24 pos=coor coor=-180,-90 breaks=30
v.mkgrid map=grid_30 box=30,30 grid=6,12 pos=coor coor=-180,-90 breaks=30

Reproject North side

  • Create new location based on EPSG #3413
# WGS 84 / NSIDC Sea Ice Polar Stereographic North
<3413> +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 \
       +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs  <>
#### NORTHERN HEMISPHERE #################
# pull in grids
v.proj in=grid_15 location=world_ll
v.proj in=grid_30 location=world_ll

g.region n=10000000 s=-10000000 w=-10000000 e=10000000 res=5000 -p

MAP=land_ocean_ice_lights

for BAND in red green blue ; do
  echo "Processing [$BAND] ..."
  r.proj -n in=$MAP.$BAND location=world_ll \
     method=cubic_f memory=300
  r.mapcalc "${MAP}_$BAND = if($MAP.$BAND > 255, 255, if($MAP.$BAND < 0, 0, int($MAP.$BAND)))"
  g.remove "$MAP.$BAND"
  g.rename "${MAP}_$BAND,$MAP.$BAND"
  r.colors $MAP.$BAND color=grey255
  echo
done

# draw
d.erase
d.rgb r=$MAP.red g=$MAP.green b=$MAP.blue
d.vect grid_15 type=boundary color=80:40:40
d.vect grid_30 type=boundary color=120:50:50

# output 
d.out.file night_on_earth_N -c size=1000,1000 format=png
d.out.file night_on_earth_N_small -c size=250,250 format=png

Reproject South side

  • Create new location based on the inverse of EPSG #3413
+proj       : stere
+datum      : wgs84
+lat_0      : -90
+lat_ts     : -70
+lon_0      : +135
+k          : 1
#### SOUTHERN HEMISPHERE #################
# pull in grids
v.proj in=grid_15 location=world_ll
v.proj in=grid_30 location=world_ll

g.region n=10000000 s=-10000000 w=-10000000 e=10000000 res=5000 -p

MAP=land_ocean_ice_lights

for BAND in red green blue ; do
  echo "Processing [$BAND] ..."
  r.proj -n in=$MAP.$BAND location=world_ll \
     method=cubic_f memory=600
  r.mapcalc "${MAP}_$BAND = if($MAP.$BAND > 255, 255, if($MAP.$BAND < 0, 0, int($MAP.$BAND)))"
  g.remove "$MAP.$BAND"
  g.rename "${MAP}_$BAND,$MAP.$BAND"
  r.colors $MAP.$BAND color=grey255
  echo
done

# remove stray NULLs
for BAND in  green blue ; do
  echo "Processing [$BAND] ..."
  g.rename "$MAP.$BAND,$MAP.$BAND.tmp1"
  r.neighbors in="$MAP.$BAND.tmp1" out="$MAP.$BAND.tmp2"
  r.patch in="$MAP.$BAND.tmp1,$MAP.$BAND.tmp2" out="$MAP.$BAND"
  g.remove "$MAP.$BAND.tmp1,$MAP.$BAND.tmp2"
  echo
done

# draw
d.erase
d.rgb r=$MAP.red g=$MAP.green b=$MAP.blue
d.vect grid_15 type=boundary color=80:40:40
d.vect grid_30 type=boundary color=120:50:50

# output
d.out.file night_on_earth_S -c size=1000,1000 format=png
d.out.file night_on_earth_S_small -c size=250,250 format=png

Result