Happy New Year! at EXPLAIN EXTENDED

Happy New Year: view of Earth from space in SQL

It's been a really long year.

The blog has not been updated during this year. I focused on improving my SQL skills and helping other folks to improve theirs.

During this year Explain Extended has become a team. We now do SQL consulting and database development. I'll write a post about it soon, but send me a message if you're impatient to wait.

Now, to the new year post.

New Year has always seemed like the most global thing to me. It does not happen to everyone at once. When I was a kid I wanted to go to space and see the New Year marching through the planet with my own eyes (and it's not that I don't want it anymore).

Blue Marble

Now we have satellite images and such, but with a little effort we can see how our planet looks from space using SQL.

To do this, we need some map data, basic math and pretty simple SQL. I'll use PostgreSQL for that.

Data

We would need some basic outlines of continents and major islands. A brief googling tells that kind folks at http://www.naturalearthdata.com/ have made it available in public domain. We would just need to go there and download the land area polygons in ESRI Shapefile format by this link.

To load that into the database I would have to create a table:

CREATE TABLE
        land
        (
        id INT PRIMARY KEY,
        path PATH NOT NULL
        );

and write a simple python script to convert ESRI files to something PostgreSQL would understand:

import shapefile

sf = shapefile.Reader("land/ne_110m_land")
print "INSERT\nINTO\tland\nVALUES\t\n" + ",\n".join(
        [
                "\t(" + str(i + 1) + ", '" + ", ".join(
                [
                        "(" + ", ".join(
                        [
                                str(round(c, 2))
                        for c in p
                        ]
                        ) +
                        ")"
                for p in b.points[:-1]
                ]
                ) +
                "'::PATH)"
        for i, b in enumerate(sf.shapes())
        ]
) + ";\n"

This gives us 127 polygons of 5016 points in total. It's quite a large query so if you can't reproduce the steps above see it on Gist here.

Showing the map

To show the map we would use some kind of a map projection. I'm not a map expert so I had to spent some time in research but it turned out that General Perspective Projection is what we need. As the Wikipedia article on the subject says, "when the Earth is photographed from space, the camera records the view as a perspective projection". That's pretty fine with me.

To render the Earth view on a "grid" we would need to provide some parameters: our coordinates, elevation, angle of view and screen resolution. Each point of the grid would correspond to some point on the Earth surface (identified by its coordinates).

If we were looking from above the North pole, the coordinates would map in a pretty straightforward way. It's simplest to work with the polar coordinates from the center of grid. Thus, φ maps directly into longitude, and R maps into the angle between the grid point and the point of view: β = ATAN(TAN(angle_of_view) * (grid_width / R)).

The point of view, the Earth center and the observed point on Earth form a triangle of which we know the two sides (distance from the center of the Earth and the Earth's radius) and the non-included angle β. We are interested in the angle between the two sides which would map into the complement to the latitude. The formula for this is γ = ASIN((height + radius) / radius * SIN(β)).

If we are not on the pole we would need to translate the points. Now, φ becomes bearing and γ becomes arc distance from our location. The translation formulae are:

LAT2 = ASIN(SIN(LAT1) * COS(γ) + COS(LAT1) * SIN(γ) * cos(φ)) LON2 = LON1 + ATAN2(SIN(φ) * SIN(γ) * COS(LAT1), COS(γ) − SIN(LAT1) * SIN(LAT2))

Now, we just need to put it into the query and observe our Earth from above Moscow:

Show the query

WITH    data (w, bmax, r, h, ilat, ilon) AS
        (
        VALUES
        (40, RADIANS(20), 6731.0, 13000.0, RADIANS(55.75), RADIANS(37.61))
        ),
        grid AS
        (
        SELECT  *
        FROM    (
                SELECT  *, generate_series(-w, w) AS x
                FROM    (
                        SELECT  *, (h + r) / r AS c,
                                generate_series(-w, w) AS y
                        FROM    data
                        ) dx
                ) dy
        ),
        points AS
        (
        SELECT  x, y, d, lat, lon
        FROM    (
                SELECT  *,
                        DEGREES(rlat) AS lat,
                        DEGREES(rlon) AS lon
                FROM    (
                        SELECT  *,
                                ilon + ATAN2(SIN(bearing) * SIN(arc) * COS(ilat), COS(arc) - SIN(ilat) * SIN(rlat)) AS rlon
                        FROM    (
                                SELECT  *,
                                        ASIN(SIN(ilat) * COS(arc) + COS(ilat) * SIN(arc) * COS(bearing)) AS rlat
                                FROM    (
                                        SELECT  *,
                                                ASIN(CASE WHEN ABS(d) < 1 THEN d END) AS arc
                                        FROM    (
                                                SELECT  *,
                                                        c * SIN(beta) AS d
                                                FROM    (
                                                        SELECT  *,
                                                                ATAN((l / w) * TAN(bmax)) AS beta
                                                        FROM    (
                                                                SELECT  *,
                                                                        ATAN2(x, -y) AS bearing,
                                                                        SQRT(X*X + Y*Y) AS l
                                                                FROM    grid
                                                                ) q
                                                        ) q
                                                ) q
                                        ) q
                                ) q
                        ) q
                ) q
        ),
        lands AS
        (
        SELECT  *,
                EXISTS
                (
                SELECT  NULL
                FROM    land
                WHERE   POINT(lon, lat) <@ path::POLYGON
                ) AS land,
                d BETWEEN 0.97 AND 1 AS border,
                x = 0 AND y = 0 AS my
        FROM    points
        )
SELECT  STRING_AGG
                (
                CASE
                WHEN my THEN
                        'x'
                WHEN border THEN
                        '*'
                WHEN land THEN
                        '#'
                ELSE ' '
                END,
                ''
                ORDER BY x
                )
FROM    lands
GROUP BY
        y
ORDER BY
        y

string_agg

*****************
**** ****
**##### # **
**########## **
*################ *
##################
*################### *
#####################
*############### # # *
######### ###### #### #####
######## ### ## #######
######## #### # # ##### ####
###### ## # # # ###### ###
######## # ###### # ########
############# # #### ### #######
* ############ ##### ### ######## *
######### ### ## ########
######## #### ######## # ##
* # ### ####### ######### ### # *
# ## ############ ################ #
# ## ############# # ################# #
* ## ############# ####################### # *
# ########### ##################### #
* ###### # # ##################### # *
* # # ###################### *
# ######################## #
* ######################### *
## # # ######################
# # ####################### #
* # # ###################### # *
* #### ########################### *
####### #################################
###### # #################################
* ##### ###################################### *
* ## #### ######################################### *
* # # # ## ######################################### *
* ## ## ###################################### *
* ## # # ######################################## #*
* # # # ######################################### *
* ################x################################### *
* #################################################### *
* ## ################################################### *
* ######################################################## *
* ######################################################## *
* #### # ############################################### *
* ### # ####### ## ##### ############################# *
### # ## ##### ###### ##############################
# ####### # #### ### # #############################
* ############## #### # #### ########################### *
* ################ # # ########### #################### ### *
############### # ########### ################### #
################# ################################ ###
* ################## ############################ #*
#################### # ############################
* ######################## ############################ ##*
* ##################################### ######### ######## *
##################################### ## ### #######
* ########################### ######### # #### *
############################ ######### #### ####
########################### ############ ###
* ############################ ############ ## *
########################## ##########
######################### ########
* ####################### ###### *
###################### ####
#########################
###################### ###
###########################
##########################
#########################
* ####################### *
####################
* ################### *
##################
*################## *
**############### **
**############ ## **
****###### ##****
*****************

81 rows fetched in 0.0009s (1.5156s)

or from above New York City:

Show the query

WITH    data (w, bmax, r, h, ilat, ilon) AS
        (
        VALUES
        (40, RADIANS(20), 6731.0, 13000.0, RADIANS(40.664167), RADIANS(-73.938611))
        ),
        grid AS
        (
        SELECT  *
        FROM    (
                SELECT  *, generate_series(-w, w) AS x
                FROM    (
                        SELECT  *, (h + r) / r AS c,
                                generate_series(-w, w) AS y
                        FROM    data
                        ) dx
                ) dy
        ),
        points AS
        (
        SELECT  x, y, d, lat, lon
        FROM    (
                SELECT  *,
                        DEGREES(rlat) AS lat,
                        DEGREES(rlon) AS lon
                FROM    (
                        SELECT  *,
                                ilon + ATAN2(SIN(bearing) * SIN(arc) * COS(ilat), COS(arc) - SIN(ilat) * SIN(rlat)) AS rlon
                        FROM    (
                                SELECT  *,
                                        ASIN(SIN(ilat) * COS(arc) + COS(ilat) * SIN(arc) * COS(bearing)) AS rlat
                                FROM    (
                                        SELECT  *,
                                                ASIN(CASE WHEN ABS(d) < 1 THEN d END) AS arc
                                        FROM    (
                                                SELECT  *,
                                                        c * SIN(beta) AS d
                                                FROM    (
                                                        SELECT  *,
                                                                ATAN((l / w) * TAN(bmax)) AS beta
                                                        FROM    (
                                                                SELECT  *,
                                                                        ATAN2(x, -y) AS bearing,
                                                                        SQRT(X*X + Y*Y) AS l
                                                                FROM    grid
                                                                ) q
                                                        ) q
                                                ) q
                                        ) q
                                ) q
                        ) q
                ) q
        ),
        lands AS
        (
        SELECT  *,
                EXISTS
                (
                SELECT  NULL
                FROM    land
                WHERE   POINT(lon, lat) <@ path::POLYGON
                ) AS land,
                d BETWEEN 0.97 AND 1 AS border,
                x = 0 AND y = 0 AS my
        FROM    points
        )
SELECT  STRING_AGG
                (
                CASE
                WHEN my THEN
                        'x'
                WHEN border THEN
                        '*'
                WHEN land THEN
                        '#'
                ELSE ' '
                END,
                ''
                ORDER BY x
                )
FROM    lands
GROUP BY
        y
ORDER BY
        y

string_agg

*****************
**** #######****
** ###########**
** ######## ####**
* ## ## #######*
# ########
* # #######*
#######
* # #######*
## ############
# ## #############
### #########
#### # # ### #####
## #### # ### ### #####
######## ### ###### #### ##### #
* # ######## ### ######## ## ##########*
# ####### ## ###### ######### ## ########
######## ## # ####### ########
* ############ ## ####### ## ##########*
######## ## ### ###### ### ## #### #
########### #### #### ##########
* ############## # ## #### ### ### #*
################# ## #### # # #### #
* ############### ## ### ## ##*
* ############## ## ## ## ## # ##*
############ ## # ## #####
* ############# ### #### ######*
############# ### # #### ######
############# ###### ### ######
* ################ ####### ##########*
* ################# ######### ########*
################## ######### ########
############################ # ########
* ######################### # ## #########*
* ######################## # ## ########*
* ############################# ########*
* ########################## # ########*
* ######################## ########*
* ######################## ########*
* ######################x ########*
* ###################### ########*
* # #################### ########*
* #################### #######*
* ################### #######*
* # ################ ###### *
* # ######## ###### ###### *
###### # ####
##### ## ###
* ##### # # *
* #### *
####
#### ### #
* #### ## # *
###### ###
* ## *
* ###### *
##
* # *
# ###### ##
#### ###########
* ################# *
##################
################## # ####
* ############################ *
############################
############################
###########################
##########################
#########################
#######################
* ####################### *
####################
* ################ *
#############
* ############ *
** ###########**
** #########**
**** #####****
*****************

81 rows fetched in 0.0010s (1.4531s)

or from above Christmas Island:

Show the query

WITH    data (w, bmax, r, h, ilat, ilon) AS
        (
        VALUES
        (40, RADIANS(20), 6731.0, 13000.0, RADIANS(-10.483333), RADIANS(105.633333))
        ),
        grid AS
        (
        SELECT  *
        FROM    (
                SELECT  *, generate_series(-w, w) AS x
                FROM    (
                        SELECT  *, (h + r) / r AS c,
                                generate_series(-w, w) AS y
                        FROM    data
                        ) dx
                ) dy
        ),
        points AS
        (
        SELECT  x, y, d, lat, lon
        FROM    (
                SELECT  *,
                        DEGREES(rlat) AS lat,
                        DEGREES(rlon) AS lon
                FROM    (
                        SELECT  *,
                                ilon + ATAN2(SIN(bearing) * SIN(arc) * COS(ilat), COS(arc) - SIN(ilat) * SIN(rlat)) AS rlon
                        FROM    (
                                SELECT  *,
                                        ASIN(SIN(ilat) * COS(arc) + COS(ilat) * SIN(arc) * COS(bearing)) AS rlat
                                FROM    (
                                        SELECT  *,
                                                ASIN(CASE WHEN ABS(d) < 1 THEN d END) AS arc
                                        FROM    (
                                                SELECT  *,
                                                        c * SIN(beta) AS d
                                                FROM    (
                                                        SELECT  *,
                                                                ATAN((l / w) * TAN(bmax)) AS beta
                                                        FROM    (
                                                                SELECT  *,
                                                                        ATAN2(x, -y) AS bearing,
                                                                        SQRT(X*X + Y*Y) AS l
                                                                FROM    grid
                                                                ) q
                                                        ) q
                                                ) q
                                        ) q
                                ) q
                        ) q
                ) q
        ),
        lands AS
        (
        SELECT  *,
                EXISTS
                (
                SELECT  NULL
                FROM    land
                WHERE   POINT(lon, lat) <@ path::POLYGON
                ) AS land,
                d BETWEEN 0.97 AND 1 AS border,
                x = 0 AND y = 0 AS my
        FROM    points
        )
SELECT  STRING_AGG
                (
                CASE
                WHEN my THEN
                        'x'
                WHEN border THEN
                        '*'
                WHEN land THEN
                        '#'
                ELSE ' '
                END,
                ''
                ORDER BY x
                )
FROM    lands
GROUP BY
        y
ORDER BY
        y

string_agg

*****************
****#############****
**#####################**
**######################### **
*############################## *
##################################
*################################## *
####################################
* ################################### ## *
############################### ###
# ############################## #
# ################################# # ##
# ################################# ####
##################################### ##
#######################################
########################################
#######################################
## ####################################
### # ############################# #
### # ########## ########## #
####### ####### ######## ##
* ####### ###### ######## # *
###### ##### ######### #
# ##### #### ####### #
# ### #### ####### #
# ### ### # ###
## # ## # #
### # # # #
## # # # #
###### # ## ###
##### # # ##
##### ## # #####
##### ### ##### ###
##### ### ##### # #
##### ### ##### ## ##
#### ## # ###
### # ####
### #### #### #
### # # ##### #
### x # ## #
## # #
### ####
### # #### #
### ####### #
#### ######### ##
#### ## ##############
### ## #################
### ## ####################
### ## ####################
## ## #################### #
# ## #####################
## # ####################
## # ####################
## ####################
# ##### #########
## #########
# #######
# #####
# ####

* *
# #

* # *
#





* *

* ############# *
#### ################
* ##### ################## *
**######################### **
**################## **
****########### ****
*****************

81 rows fetched in 0.0009s (1.4531s)

Happy New Year!

Previous New Year posts:

Written by Quassnoi

December 31st, 2012 at 11:00 pm

Posted in Uncategorized

Tagged with ASCII art, map projection, PostgreSQL

2 Responses to 'Happy New Year: view of Earth from space in SQL'

Subscribe to comments with RSS

Nice query! What tool do you use to show the output of your queries? I installed PostgreSQL on Win 7, and it comes with psql which is a basic command line tool to run queries.

Tyler

17 Jan 13 at 08:40
@Tyler: that’s what I use, psql. I just copy its output and paste into the post.

Quassnoi

17 Jan 13 at 10:25

string_agg

*****************
**** ****
**##### # **
**########## **
*################ *
##################
*################### *
#####################
*############### # # *
######### ###### #### #####
######## ### ## #######
######## #### # # ##### ####
###### ## # # # ###### ###
######## # ###### # ########
############# # #### ### #######
* ############ ##### ### ######## *
######### ### ## ########
######## #### ######## # ##
* # ### ####### ######### ### # *
# ## ############ ################ #
# ## ############# # ################# #
* ## ############# ####################### # *
# ########### ##################### #
* ###### # # ##################### # *
* # # ###################### *
# ######################## #
* ######################### *
## # # ######################
# # ####################### #
* # # ###################### # *
* #### ########################### *
####### #################################
###### # #################################
* ##### ###################################### *
* ## #### ######################################### *
* # # # ## ######################################### *
* ## ## ###################################### *
* ## # # ######################################## #*
* # # # ######################################### *
* ################x################################### *
* #################################################### *
* ## ################################################### *
* ######################################################## *
* ######################################################## *
* #### # ############################################### *
* ### # ####### ## ##### ############################# *
### # ## ##### ###### ##############################
# ####### # #### ### # #############################
* ############## #### # #### ########################### *
* ################ # # ########### #################### ### *
############### # ########### ################### #
################# ################################ ###
* ################## ############################ #*
#################### # ############################
* ######################## ############################ ##*
* ##################################### ######### ######## *
##################################### ## ### #######
* ########################### ######### # #### *
############################ ######### #### ####
########################### ############ ###
* ############################ ############ ## *
########################## ##########
######################### ########
* ####################### ###### *
###################### ####
#########################
###################### ###
###########################
##########################
#########################
* ####################### *
####################
* ################### *
##################
*################## *
**############### **
**############ ## **
****###### ##****
*****************

81 rows fetched in 0.0009s (1.5156s)

string_agg

*****************
**** #######****
** ###########**
** ######## ####**
* ## ## #######*
# ########
* # #######*
#######
* # #######*
## ############
# ## #############
### #########
#### # # ### #####
## #### # ### ### #####
######## ### ###### #### ##### #
* # ######## ### ######## ## ##########*
# ####### ## ###### ######### ## ########
######## ## # ####### ########
* ############ ## ####### ## ##########*
######## ## ### ###### ### ## #### #
########### #### #### ##########
* ############## # ## #### ### ### #*
################# ## #### # # #### #
* ############### ## ### ## ##*
* ############## ## ## ## ## # ##*
############ ## # ## #####
* ############# ### #### ######*
############# ### # #### ######
############# ###### ### ######
* ################ ####### ##########*
* ################# ######### ########*
################## ######### ########
############################ # ########
* ######################### # ## #########*
* ######################## # ## ########*
* ############################# ########*
* ########################## # ########*
* ######################## ########*
* ######################## ########*
* ######################x ########*
* ###################### ########*
* # #################### ########*
* #################### #######*
* ################### #######*
* # ################ ###### *
* # ######## ###### ###### *
###### # ####
##### ## ###
* ##### # # *
* #### *
####
#### ### #
* #### ## # *
###### ###
* ## *
* ###### *
##
* # *
# ###### ##
#### ###########
* ################# *
##################
################## # ####
* ############################ *
############################
############################
###########################
##########################
#########################
#######################
* ####################### *
####################
* ################ *
#############
* ############ *
** ###########**
** #########**
**** #####****
*****************

81 rows fetched in 0.0010s (1.4531s)

string_agg

*****************
****#############****
**#####################**
**######################### **
*############################## *
##################################
*################################## *
####################################
* ################################### ## *
############################### ###
# ############################## #
# ################################# # ##
# ################################# ####
##################################### ##
#######################################
########################################
#######################################
## ####################################
### # ############################# #
### # ########## ########## #
####### ####### ######## ##
* ####### ###### ######## # *
###### ##### ######### #
# ##### #### ####### #
# ### #### ####### #
# ### ### # ###
## # ## # #
### # # # #
## # # # #
###### # ## ###
##### # # ##
##### ## # #####
##### ### ##### ###
##### ### ##### # #
##### ### ##### ## ##
#### ## # ###
### # ####
### #### #### #
### # # ##### #
### x # ## #
## # #
### ####
### # #### #
### ####### #
#### ######### ##
#### ## ##############
### ## #################
### ## ####################
### ## ####################
## ## #################### #
# ## #####################
## # ####################
## # ####################
## ####################
# ##### #########
## #########
# #######
# #####
# ####

* *
# #

* # *
#





* *

* ############# *
#### ################
* ##### ################## *
**######################### **
**################## **
****########### ****
*****************

81 rows fetched in 0.0009s (1.4531s)

EXPLAIN EXTENDED

Happy New Year: view of Earth from space in SQL

Data

Showing the map

2 Responses to 'Happy New Year: view of Earth from space in SQL'

Leave a Reply

Subscribe

Subscribe by email

Contacts

Should I?

Recent articles

Calendar

Archives

Categories

Stack Overflow

EXPLAIN EXTENDED

Happy New Year: view of Earth from space in SQL

Data

Showing the map

Share this:

2 Responses to 'Happy New Year: view of Earth from space in SQL'

Leave a Reply

Subscribe

Subscribe by email

Contacts

Should I?

Recent articles

Calendar

Archives

Categories

Stack Overflow