From Lat/Lon to UTM zone in Spatialite

There are several PostGIS functions floating around to calculate the UTM zone EPSG code for points in Latitude/Longitude WGS84. However, Spatialite, based on Sqlite, does not support user created functions. So how can we get the same results in a Spatialite database?

[Correction: Spatialite does, of course, have built in functions, and users can create their own functions using the C API. However - unlike PostGIS - there is no support for creating SQL (or other procedural languages) based functions. SQlite does not have the equivalent of the CREATE FUNCTION command in SQL]

The EPSG codes for UTM zones can be calculated as (Longitude+186)/6 + 32600. In the southern hemisphere, it’s (Long+186)/6 + 32700. The last two digits of each EPSG code are the UTM zone. The third digit, 6 or 7 indicates N or S. We can use the Sqlite CASE WHEN…ELSE…END construct to get the correct EPSG code for a table of point locations as follows:

We begin with the freely available geonames cities15000 file. This is a tab delimited text file of some 22,000 cities (with population over 15,000) around the world. I linked to the table, (making a virtual table in Spatialite), and used the data to create a real table with more convenient column names. I then made the table spatial, and used the latitude and longitude columns to set the Geometry for all cities.

id integer primary key autoincrement, geoname_code integer, Name text, Name_en text,
lat double, lon double, Country_ISO text, pop integer, Country text);


(geoname_code, Name, Name_en, lat, lon, Country_ISO, pop, Country)


col001, col002, col003, col005, col006, col009, col015, col018

FROM cities15000;

Now setup the Geometry column

SELECT AddGeometryColumn('cities','Geometry',4326,'POINT','XY');
UPDATE cities SET Geometry=MakePoint(lon, lat, 4326);

And here’s how to construct a query to get the UTM zone for each city:

SELECT Name_en AS "City", Country AS "Continent/Country",

CASE WHEN ST_X(Geometry)=180 THEN
(CAST(ST_X(Geometry) AS INTEGER)+186)/6+32600


CASE WHEN ST_X(Geometry)=180 THEN
(CAST(ST_X(Geometry) AS INTEGER)+186)/6+32700

FROM cities
WHERE SRID(Geometry)=4326;

A few comments:

  • The final WHERE clause insures that no EPSG codes will be returned if the point layer is not in Lat/Lon WGS84 to begin with.
  • The outer CASE statement separates Northern and Southern hemispheres.
  • The inner CASE catches the edges where longitude values are 180 degrees.

Here are some results:

City          Continent/Country  UTM SRID  
------------  -----------------  ----------
Sao Jeronimo  America/Sao_Paulo  32722     
Jeremoabo     America/Bahia      32724     
Saint-Jerome  America/Montreal   32618     
Aguada de Pa  America/Havana     32617     
Esbjerg       Europe/Copenhagen  32632     
Jerez de la   Europe/Madrid      32630     
Jeremie       America/Port-au-P  32618     
Jerusalem     Asia/Jerusalem     32636     
West Jerusal  Asia/Jerusalem     32636     
Manjeri       Asia/Kolkata       32643     
Borujerd      Asia/Tehran        32639     
Jerada        Africa/Casablanca  32630     
San Jeronimo  America/Mexico_Ci  32615     
Jerez de Gar  America/Mexico_Ci  32613     
Jerantut      Asia/Kuala_Lumpur  32648     
Ijero-Ekiti   Africa/Lagos       32631     
Oud-Beijerla  Europe/Amsterdam   32631     
Steinkjer     Europe/Oslo        32632     
Jericho       Asia/Hebron        32636     
East Jerusal  Asia/Hebron        32636     
Oum Hadjer    Africa/Ndjamena    32634     
Tajerouine    Africa/Tunis       32632     
Njeru         Africa/Kampala     32636     
Jersey City   America/New_York   32618     
Las Tejerias  America/Caracas    32619

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>