source: DPPP/ukcip02_server/trunk/bng2wgs84.py @ 3360

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/DPPP/ukcip02_server/trunk/bng2wgs84.py@5243
Revision 3360, 1.7 KB checked in by spascoe, 12 years ago (diff)

Quick import of code I have been working on to server ukcip02 data.
This code may well move soon.

Line 
1#!/usr/bin/env/python2.4
2"""
3Proof of concept conversion of OSGB36 grid to WGS84
4
5@author: Stephen Pascoe
6
7"""
8
9import sys, os
10
11from osgeo import gdal, osr
12from osgeo.gdalconst import *
13import cdms
14
15import numpy
16from numpy.core import ma
17
18import matplotlib
19matplotlib.rcParams['numerix'] = 'numpy'
20from matplotlib.cm import get_cmap
21from matplotlib.colors import normalize
22
23
24
25dataFile = sys.argv[1]
26f = cdms.open(dataFile)
27
28# Get first time step, first level
29varName = f.listvariables()[0]
30var = f[varName][0]
31
32a = ma.array(var.getValue(), mask=var.mask())
33
34# Get min and max
35a_min = numpy.min(a)
36a_max = numpy.max(a)
37
38# Convert to RGBA array
39cmap = get_cmap()
40a_colour = cmap(normalize(a_min, a_max)(a))
41
42# Create tif of same dimensions
43y, x = a.shape
44d = gdal.GetDriverByName('GTIFF').Create('out.tif', x, y, 4, GDT_Byte)
45
46for i in range(4):
47    d.GetRasterBand(i+1).WriteArray(a_colour[:,:,i]*255)
48
49
50srs = osr.SpatialReference()
51srs.ImportFromEPSG(27700)
52
53print '--'
54print srs.ExportToWkt()
55print srs.ExportToProj4()
56print '--'
57
58northings, eastings = var.getAxisList()
59# Have I got these the right way round?
60assert(northings.id == 'northings')
61
62n_min = min(northings[0], northings[-1])
63n_max = max(northings[0], northings[-1])
64e_min = min(eastings[0], eastings[-1])
65e_max = max(eastings[0], eastings[-1])
66
67d.SetGCPs([gdal.GCP(eastings[0], northings[-1], 0, 0.5,0.5),
68           gdal.GCP(eastings[-1], northings[-1], 0, x+0.5,0.5),
69           gdal.GCP(eastings[0], northings[0], 0, 0.5,y+0.5),
70           gdal.GCP(eastings[-1], northings[0], 0, x+0.5,y+0.5)],
71          srs.ExportToWkt())
72del d
73
74# Warping bindings are incomplete therefore use os.system
75os.system('gdalwarp -q -t_srs EPSG:4326 -te -15 45 5 60 out.tif out2.tif')
76
Note: See TracBrowser for help on using the repository browser.