Earth at night: Difference between revisions

From GRASS-Wiki
Jump to navigation Jump to search
(+img)
m (→‎Import: fixed map name according to link fix)
 
(16 intermediate revisions by 3 users not shown)
Line 1: Line 1:
== Import NASA's Night on Earth and Reproject in Polar Sterographic ==
== Import NASA's Night on Earth and Reproject in Polar Stereographic ==


''This micro-tutorial requires GRASS 6.5''
''This micro-tutorial requires GRASS 6.5''
Line 6: Line 6:


Download
Download
http://veimages.gsfc.nasa.gov//1438/land_ocean_ice_lights_8192.tif
http://eoimages.gsfc.nasa.gov/ve/1438/land_ocean_ice_lights_2048.tif
from
from NASA's
http://visibleearth.nasa.gov/view_detail.php?id=1438
[http://visibleearth.nasa.gov/view_set.php?categoryID=2363 Visible Earth] website.


=== Import ===
=== Import ===
  # create simple XY location
* Create a WGS84 Lat/Lon location (EPSG # 4326)
  r.in.gdal in=land_ocean_ice_lights_8192.tif out=land_ocean_ice_lights
 
  # 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_2048.tif out=land_ocean_ice_lights
   
   
  for BAND in red green blue ; do
  for BAND in red green blue ; do
Line 18: Line 24:
  done
  done
   
   
  #FIXME
  g.region n=90 s=45 w=-180 e=180
#hack PROJ WIND and cellhd to be 3 not 0
#run g.setproj to set datum == wgs84 and create PROJ_INFO etc.
# now this location is lat/lon
   
   
  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_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
  v.mkgrid map=grid_30 box=30,30 grid=6,12 pos=coor coor=-180,-90 breaks=30


=== Reproject North side ===
=== 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 #################
  #### NORTHERN HEMISPHERE #################
#create new loc'n 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  <>
  # pull in grids
  # pull in grids
  v.proj in=grid_15 location=xy_to_ll
  v.proj in=grid_15 location=world_ll
  v.proj in=grid_30 location=xy_to_ll
  v.proj in=grid_30 location=world_ll
   
   
d.grid -w 30 color=red
  g.region n=10000000 s=-10000000 w=-10000000 e=10000000 res=5000 -p
  g.region n=10000000 s=-10000000 w=-10000000 e=10000000 res=5000 -p
   
   
Line 45: Line 48:
  for BAND in red green blue ; do
  for BAND in red green blue ; do
   echo "Processing [$BAND] ..."
   echo "Processing [$BAND] ..."
   r.proj -n in=$MAP.$BAND location=xy_to_ll \
   r.proj -n in=$MAP.$BAND location=world_ll \
       method=cubic_f memory=600
       method=cubic_f memory=300
   r.mapcalc "${MAP}_$BAND = if($MAP.$BAND > 255, 255, if($MAP.$BAND < 0, 0, int($MAP.$BAND)))"
   r.mapcalc "${MAP}_$BAND = if($MAP.$BAND > 255, 255, if($MAP.$BAND < 0, 0, int($MAP.$BAND)))"
   g.remove "$MAP.$BAND"
   g.remove "$MAP.$BAND"
Line 54: Line 57:
  done
  done
   
   
# draw
  d.erase
  d.erase
  d.rgb r=$MAP.red g=$MAP.green b=$MAP.blue
  d.rgb r=$MAP.red g=$MAP.green b=$MAP.blue
Line 59: Line 63:
  d.vect grid_30 type=boundary color=120:50:50
  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 -c size=1000,1000 format=png
  d.out.file night_on_earth_N_small -c size=250,250 format=png
  d.out.file night_on_earth_N_small -c size=250,250 format=png


=== Reproject South side ===
=== 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 #################
  #### SOUTHERN HEMISPHERE #################
  #create new loc'n based on inverse of epsg #3413
  # pull in grids
  #proj       : stere
  v.proj in=grid_15 location=world_ll
  #datum      : wgs84
  v.proj in=grid_30 location=world_ll
#lat_0      : 90
#lat_ts    : 70
#lon_0      : +135
#k          : 1
   
   
  v.proj in=grid_15 location=xy_to_ll
  g.region n=10000000 s=-10000000 w=-10000000 e=10000000 res=5000 -p
v.proj in=grid_30 location=xy_to_ll
   
   
g.region n=10000000 s=-10000000 w=-10000000 e=10000000 res=5000 -p
  MAP=land_ocean_ice_lights
  MAP=land_ocean_ice_lights
   
   
  for BAND in red green blue ; do
  for BAND in red green blue ; do
   echo "Processing [$BAND] ..."
   echo "Processing [$BAND] ..."
   r.proj -n in=$MAP.$BAND location=xy_to_ll \
   r.proj -n in=$MAP.$BAND location=world_ll \
       method=cubic_f memory=600
       method=cubic_f memory=600
   r.mapcalc "${MAP}_$BAND = if($MAP.$BAND > 255, 255, if($MAP.$BAND < 0, 0, int($MAP.$BAND)))"
   r.mapcalc "${MAP}_$BAND = if($MAP.$BAND > 255, 255, if($MAP.$BAND < 0, 0, int($MAP.$BAND)))"
Line 100: Line 107:
  done
  done
   
   
# draw
  d.erase
  d.erase
  d.rgb r=$MAP.red g=$MAP.green b=$MAP.blue
  d.rgb r=$MAP.red g=$MAP.green b=$MAP.blue
Line 105: Line 113:
  d.vect grid_30 type=boundary color=120:50:50
  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 -c size=1000,1000 format=png
  d.out.file night_on_earth_S_small -c size=250,250 format=png
  d.out.file night_on_earth_S_small -c size=250,250 format=png


=== Result ===
=== Result ===

Latest revision as of 12:51, 15 September 2011

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_2048.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

Retrieved from "https://grasswiki.osgeo.org/w/index.php?title=Earth_at_night&oldid=14044"