Earth at night

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

This micro-tutorial requires GRASS 6.5

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=xy_to_ll
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

for BAND in red green blue ; do
echo "Processing [\$BAND] ..."
r.proj -n in=\$MAP.\$BAND location=xy_to_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

```#### SOUTHERN HEMISPHERE #################
#create new loc'n based on inverse of epsg #3413
#proj       : stere
#datum      : wgs84
#lat_0      : 90
#lat_ts     : 70
#lon_0      : +135
#k          : 1

# pull in grids
v.proj in=grid_15 location=xy_to_ll
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

for BAND in red green blue ; do
echo "Processing [\$BAND] ..."
r.proj -n in=\$MAP.\$BAND location=xy_to_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
```