Classifying point coordinates to countries or continents?

Question

I have about 50 000 points (latitude and longitude), and I would like to find out in which country each of the points is located. It does not need to be super precise, though.

(Or actually I only need to classify them to larger geographical areas (North America, South America, Europe, Africa, Middle East, Ex-USSR, Asia, Oceania), but I assume classifying the points to countries and then the countries to larger geographical area is the way to do it.)

How should I proceed?

Preferably the tool to do this should be free, run on Linux, and be a command-line tool, or e.g. Python or R program (I can write the program if I find some starting points or hints) or some library that I could use. I have no experience in GIS software, and I would prefer a programmable or a command line solution.

Answer 1

Using SAGA-GIS you can use the tool Add Polygon Attributes to Points

On the commandline (linux or windows):

$ saga_cmd shapes_points 10

library path:   C:\Program Files\SAGA-GIS\Modules
library name:   shapes_points
module name :   Add Polygon Attributes to Points
author      :   O.Conrad (c) 2009

Usage: 10 -INPUT <str> [-OUTPUT <str>] -POLYGONS <str> [-FIELD <str>]
  -INPUT:<str>          Points
        Shapes (input)
  -OUTPUT:<str>         Result
        Shapes (optional output)
  -POLYGONS:<str>       Polygons
        Shapes (input)
  -FIELD:<str>          Attribute
        Table field

Answer 2

In R, you can use my geonames package to access the www.geonames.org API:

> library(geonames)
Loading required package: rjson
Warning message:
In fun(libname, pkgname) :
  No geonamesUsername set. See http://geonames.wordpress.com/2010/03/16/ddos-part-ii/ and set one with options(geonamesUsername="foo") for some services to work
> GNcountryCode(lat=20,lng=0)$countryCode
[1] "ML"
> GNcountryCode(lat=30,lng=0)$countryCode
[1] "DZ"
> GNcountryCode(lat=40,lng=0)$countryCode
[1] "ES"
> GNcountryCode(lat=50,lng=0)$countryCode
Error in getJson("countryCode", list(lat = lat, lng = lng, radius = radius,  : 
  error code 15 from server: no country code found

I think 50,0 is in the sea. You'd have to trap this error in a loop. The function is not vectorised so you can't feed it all 50,000 points without a loop anyway.

geonames is on CRAN, so install.packages("geonames") will get it for you.

Answer 3

Following warrieka's hint, I managed the following:

# my_pointfile.txt looks the following, with tabs as separators:
# lon      lat
# 176      63
# -85      72.5
# etc.
#
# Convert text table to Saga GIS table
# (2 - Import Text Table with Numbers only)
saga_cmd libio_table 2 -TABLE my_table -SEPARATOR 0 -HEADLINE -FILENAME my_pointfile.txt

# Convert Saga table to pointfile
# (0 - Convert Table to Points)
saga_cmd libshapes_points 0 -TABLE my_table -Y 1 -POINTS my_points
# (I don't understand why -X 1 -Y 2 does not work. -Y 1 seems to work,
# but kinda funny to say -Y 1 to mean that x-coordinates are taken from
# column 1, and y from column 2.)

# Add labels from polygons to points
#   Country data is from:
#   http://www.naturalearthdata.com/downloads/50m-cultural-vectors/
# (10 - Add Polygon Attributes to Points)
# subcontinent:
saga_cmd libshapes_points 10 -INPUT my_points.shp -OUTPUT subregion -POLYGONS ne_50m_admin_0_countries.shp -FIELD subregion
# country or territory:
saga_cmd libshapes_points 10 -INPUT my_points.shp -OUTPUT country -POLYGONS ne_50m_admin_0_countries.shp -FIELD admin

# Convert shapefile (which contains points and their labels) to text
ogr2ogr -f CSV subregion.csv subregion.shp
ogr2ogr -f CSV country.csv country.shp

I didn't manage to put the subcontinent labels and country name labels in the same shapefile or .csv file. And, some points get an empty string as their label, as they are not inside any country borders (although they might be very close). I didn't find a feature to get the label from the nearest polygon.