Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question
add comment

3 Answers

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
share|improve this answer
1  
Get the polygons with geographical area's from naturalearthdata.com if you dont have them. –  warrieka Jan 29 '13 at 13:28
add comment

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.

share|improve this answer
    
Yes, that would be nice and simple, but looking at the pricing page of geonames.org, even for paying customers they seem to serve only some thousand requests (credits) per hour, so I'd guess they would not like to give me 50 000 queries ( geonames.org/commercial-webservices.html ). –  Sampo Smolander Jan 29 '13 at 15:06
1  
You can spread the geocoding process over a few days. We have geocoding processes and set a limit on the number per day using yahoo,google,mapquest etc. –  Mapperz Jan 29 '13 at 15:38
add comment
up vote 1 down vote accepted

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.