PostGIS How-to: NetCDF to PostGIS
Posted on 15 November 2011
Disclaimer: I am not sure if this is the quickest way of importing NetCDF files to PostGIS, but it works, and for me that is the most important. If you have any other suggestions for improvements, please post them in the comment field!
In my work we need to combine a lot of spatial data from diverse sources. This includes socio-economic, physical and climatic data. Our latest requirement is to be able to import NetCDF files containing climate model data into our PostGIS database. NetCDF (Network Common Data Format) is a multidimensional spatio-temporal data structure with the usual x, y spatial coordinates and a z coordinate or value in addition to the t variable which represents time since some start date.
Our approach is to read the NetCDF in Python using the Scipy.io NetCDF module. This makes it possible to read the NetCDF file into Python and then output it to PostGIS using the psycopg2 module.
I will paste my python script below and comment on the various parts. Please feel free to replicate this script, but please site the project paper:
Tollefsen, Andreas Forø; Håvard Strand & Halvard Buhaug (2012) PRIO-GRID: A unified spatial data structure. Journal of Peace Research 49(2): 363-374.
Script:
# NetCDF to PostGreSQL database
# Delaware precipitation and temperature data. From NetCDF to database table
# Requires Python2.6, Postgresql, Psycopg2, Scipy
# Tested using Vista 64bit.
# Import modules
import psycopg2, time, datetime
from scipy.io import netcdf
# Establish connection
db1 = psycopg2.connect("host=xxxx dbname=priogrid user=xxxx password=xxxx")
cur = db1.cursor()
# Create Table in postgis
print str(time.ctime())+ " Creating precipud19002008 table."
cur.execute("DROP TABLE IF EXISTS precipud19002008;")
cur.execute("CREATE TABLE precipud19002008 (gid serial PRIMARY KEY not null, year int, month int, lon decimal, lat decimal, prec decimal);")
# Read netcdf file using the netcdf_file class. 'r' for read.
f = netcdf.netcdf_file('/mnt/pc258/prio_grid/source/ClimateData/precipud19002008.nc', 'r')
# Create lathash. Read the f.variables['lat'].data.tolist() creates a list of lat coodinates.
print str(time.ctime())+ " Looping through lat coords."
temp = f.variables['lat'].data.tolist()
# Create list of lat coords
lathash = {}
for entry in temp:
lathash[temp.index(entry)] = entry
# Create list of lon coords
print str(time.ctime())+ " Looping through long coords."
temp = f.variables['lon'].data.tolist()
lonhash = {}
for entry in temp:
lonhash[temp.index(entry)] = entry
# Loop through every observation. Set timedimension and lat and long observations.
for _month in xrange(1308):
print str(time.ctime())+ " Extracting NetCDF information for month number "+str(_month+1)+" of 1308"
for _lon in xrange(720):
for _lat in xrange(360):
thisyear = int((_month)/12+1900)
thismonth = ((_month) % 12)+1
thisdate = datetime.date(thisyear,thismonth, 1)
data = [int(thisyear), int(thismonth), lonhash[_lon]-180.00, lathash[_lat], f.variables[('data')].data[_month, _lat, _lon]+32767]
cur.execute("INSERT INTO precipud19002008 (year, month, lon, lat, prec) VALUES "+str(tuple(data))+";")
print str(time.ctime())+ " Done!"
db1.commit()
cur.close()
db1.close()
So, what this scripts does is to read the NetCDF file into Python. Then loops through the month, lat, lon, creates a list of the respective values for each lat and long and insert these into the PostGIS table. This NetCDF file includes 1308 months of data. With 259200 observations per year. It takes approximately 30 seconds per month.
-
MikeT
-
Anonymous
-
Me
-
Charlie Sharpsteen