I have a 500 x 500 m resolution world map of clumping index data (He et al., 2012) and I am trying to do two things at the same time: 1) to isolate the variables per plant functional type following the GLC2000 land cover map with 1 x 1 km resolution, and 2) to rescale the new maps to 0.5 x 0.5 degree resolution.
But the maps are written in different formats:
1) gdalinfo MCI2_red_global_sza0_mask.img
Driver: ENVI/ENVI .hdr Labelled
Files: MCI2_red_global_sza0_mask.img
MCI2_red_global_sza0_mask.img.hdr
Size is 86400, 43200
Coordinate System is:
PROJCS["WGS_1984_Sinusoidal",
GEOGCS["GCS_sphere_modis",
DATUM["custom",
SPHEROID["custom",6371007.181,0.0]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Sinusoidal"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["longitude_of_center",0.0],
UNIT["Meter",1]]
Origin = (-20015109.353999998420477,10007554.676999999210238)
Pixel Size = (463.312716529999989,-463.312716529999989)
Metadata:
Band_1=Mask (Reprojection (Layer_1:glc2000_v1_1.img):MCI2_red_global_sza0.img)
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-20015109.354,10007554.677) (124d17'29.21"W, 90d 0' 0.00"N)
Lower Left (-20015109.354,-10007554.677) (154d42'44.97"E, 90d 0' 0.00"S)
Upper Right (20015109.354,10007554.677) (136d31'56.59"E, 90d 0' 0.00"N)
Lower Right (20015109.354,-10007554.677) (141d 0'27.35"W, 90d 0' 0.00"S)
Center ( 0.0000960, -0.0000480) ( 0d 0' 0.00"E, 0d 0' 0.00"S)
Band 1 Block=86400x1 Type=Byte, ColorInterp=Undefined
Description = Mask (Reprojection (Layer_1:glc2000_v1_1.img):MCI2_red_global_sza0.img)
2) gdalinfo glc2000_v1_1.img
Driver: HFA/Erdas Imagine Images (.img)
Files: glc2000_v1_1.img
glc2000_v1_1.igw
glc2000_v1_1.rrd
Size is 40320, 16353
Coordinate System is:
GEOGCS["BAD DATUM",
DATUM["BAD DATUM",
SPHEROID["WGS84",6378137,298.2570248822722]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]]
Origin = (-180.004464285700010,89.995535666300043)
Pixel Size = (0.008928571400000,-0.008928571400000)
Corner Coordinates:
Upper Left (-180.0044643, 89.9955357) (180d 0'16.07"W, 89d59'43.93"N)
Lower Left (-180.0044643, -56.0133924) (180d 0'16.07"W, 56d 0'48.21"S)
Upper Right ( 179.9955346, 89.9955357) (179d59'43.92"E, 89d59'43.93"N)
Lower Right ( 179.9955346, -56.0133924) (179d59'43.92"E, 56d 0'48.21"S)
Center ( -0.0044649, 16.9910716) ( 0d 0'16.07"W, 16d59'27.86"N)
Band 1 Block=256x4 Type=Byte, ColorInterp=Palette
Description = Layer_1
Min=1.000 Max=22.000
Minimum=1.000, Maximum=22.000, Mean=17.701, StdDev=4.927
Overviews: 10078x4086, 5038x2042, 2518x1020, 1258x509, 628x253, 313x125, 155x61
So, here it's my code:
from osgeo import gdal
import numpy as np
import sys
from array import array
import osgeo.osr as osr
# Enable exceptions
gdal.UseExceptions()
MODIS_Wkt='PROJCS["unnamed",GEOGCS["Unknown datum based upon the custom spheroid",DATUM["Not specified (based on custom spheroid)",SPHEROID["Custom spheroid",6371007.181,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Sinusoidal"],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]'
"""
A class to extract and process regions from He and Chen's Global Clumping Index data
T Quaife, 2013
"""
################################# GLC Global Class (according to LCCS terminology) ##########################################
##( 01 ) Tree Cover, broadleaved, evergreen -
## LCCS >15% tree cover, tree height >3m
## (Examples of sub-classes at regional level* :
## closed > 40% tree cove; open 15-40% tree cover)
##( 02 )Tree Cover, broadleaved, deciduous, closed
##( 03 )Tree Cover, broadleaved, deciduous, open -
## (open 15-40% tree cover)
##( 04 )Tree Cover, needle-leaved, evergreen
##( 05 )Tree Cover, needle-leaved, deciduous
##( 06 )Tree Cover, mixed leaf type
##( 07 )Tree Cover, regularly flooded, fresh water (& brackish)
##( 08 )Tree Cover, regularly flooded, saline water, -
## (daily variation of water level)
##( 09 )Mosaic: -
## Tree cover / Other natural vegetation
##( 10 )Tree Cover, burnt
##( 11 )Shrub Cover, closed-open, evergreen
## (Examples of sub-classes at reg. level *: (i) sparse tree layer)
##( 12 )Shrub Cover, closed-open, deciduous
## (Examples of sub-classes at reg. level *: (i) sparse tree layer)
##( 13 )Herbaceous Cover, closed-open
## (Examples of sub-classes at regional level *:
## (i) natural, (ii) pasture, (iii) sparse trees or shrubs)
##( 14 )Sparse Herbaceous or sparse Shrub Cover
##( 15 )Regularly flooded Shrub and/or Herbaceous Cover
##( 16 )Cultivated and managed areas
## (Examples of sub-classes at reg. level *:
## (i) terrestrial; (ii) aquatic (=flooded during cultivation), and under terrestrial:
## (iii) tree crop & shrubs (perennial),
## (iv) herbaceous crops (annual), non-irrigated, (v) herbaceous crops (annual), irrigated)
##( 17 )Mosaic: -
## Cropland / Tree Cover / Other natural vegetation
##( 18 )Mosaic: -
## Cropland / Shrub or Grass Cover
##( 19 )Bare Areas
##( 20 )Water Bodies (natural & artificial)
##( 21 )Snow and Ice (natural & artificial)
##( 22 )Artificial surfaces and associated areas
##( 23 )No data
TRANS=[]
TRANS.append(5) #lc0 - NaN
TRANS.append(0) #lc1
TRANS.append(0) #lc2
TRANS.append(0) #lc3
TRANS.append(1) #lc4
TRANS.append(1) #lc5
TRANS.append(0) #lc6
TRANS.append(0) #lc7
TRANS.append(0) #lc8
TRANS.append(0) #lc9
TRANS.append(0) #lc10
TRANS.append(4) #lc11
TRANS.append(4) #lc12
TRANS.append(2) #lc13
TRANS.append(2) #lc14
TRANS.append(2) #lc15
TRANS.append(3) #lc16
TRANS.append(3) #lc17
TRANS.append(3) #lc18
TRANS.append(5) #lc19
TRANS.append(5) #lc20
TRANS.append(5) #lc21
TRANS.append(5) #lc22
TRANS.append(5) #lc23
class globClump:
def __init__( self, filename ):
#basic GDAL set up:
self.filename=filename
self.dataSet=gdal.Open( filename, gdal.GA_ReadOnly )
self.data=self.dataSet.GetRasterBand(1)
self.dataArray=self.data.ReadAsArray()
#geographic set up:
self.wkt=self.dataSet.GetProjection()
self.proj=osr.SpatialReference( )
#n.b. the following appears not to work!
#so instead import the MODIS wkt...
#self.proj.ImportFromWkt( self.wkt )
self.proj.ImportFromWkt( MODIS_Wkt )
self.geoLL=self.proj.CloneGeogCS( )
self.transformer=osr.CoordinateTransformation( self.proj, self.geoLL )
#coordinates from GLC2000
self.glcLonMin = -180.004464285700010
self.glcLatMin = 89.995535666300043
self.glcLatMax = -56.0133924
self.glcPixelWidth = 0.0089285714 #in degress
self.glcPixelWidth = 0.0089285714 #in degress
GLCfilename = "/glusterfs/phd/users/mn811042/glc2000/Img/glc2000_v1_1.img"
self.glcDataSet=gdal.Open( GLCfilename, gdal.GA_ReadOnly )
self.glcData=self.glcDataSet.GetRasterBand(1)
self.glcDataArray=self.glcData.ReadAsArray()
transform = self.glcDataSet.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]
print 'xOrigin, yOrigin, pixelWidth, pixelHeight'
print xOrigin, yOrigin, pixelWidth, pixelHeight
sys.exit()
def getLC_fromGLC2000(self, lat, lon):
x = int( (lon-self.glcLonMin)/self.glcPixelWidth )
y = int( (self.glcLatMin-lat)/self.glcPixelWidth )
#read the data point
lcType = self.glcDataArray[y,x]
#print 'x,y,lon,lat,lcType'
#print x,y,lon,lat,lcType
return lcType
def getPointLatLon( self, row, col ):
"""
Get the lat and lon of a pixel.
Starts at 0,0.
x0/y0 = top left corrner in meters
xs/ys = x and y pixel sizes
Note xs and ys 2nd and 6th in the argument list.
"""
(x0,xs,xz,y0,yz,ys)=self.dataSet.GetGeoTransform()
#pixel coordinates in meters
pixX = col*xs + x0 + xs/2.
pixY = y0 + row*ys + ys/2.
#print (x0,xs,xz,y0,yz,ys), pixX, pixY
return self.transformer.TransformPoint( pixX, pixY )
def getLatLonAsciiMap( self, res=0.5 ):
(Y,X)=np.shape(self.dataArray)
nPft=6
outMap=np.zeros([nPft,int(180/float(res)),int(360/float(res))])
outMapN=np.zeros([nPft,int(180/float(res)),int(360/float(res))])
# 43200 86400
#print 'Y=',Y,'X=',X
#ybox= 43200
#xbox= 86400
ybox = 21600
#xbox = 50000
xbox = 28287
ybox_var= 2160
xbox_var= 2828
#for x in xrange(X):
#for x in [43200]:
for x in xrange(xbox-xbox_var,xbox+xbox_var):
print >> sys.stderr, x,
#for y in xrange(Y):
for y in xrange(ybox-ybox_var,ybox+ybox_var):
t=g.getPointLatLon( y, x )
#print x,y,t[0],t[1],self.dataArray[y,x]
data=self.dataArray[y,x]
if t[1] > self.glcLatMax:
lc = self.getLC_fromGLC2000(t[1],t[0])
pft = TRANS[lc]
else:
lc = 23
pft = TRANS[lc]
#print 'x,y,lat,lon,ci,lc,pft'
#print x,y,t[1],t[0],self.dataArray[y,x],lc,pft
if data > 0.:
yIndx=int(np.floor((t[1]+90.)/float(res)))
xIndx=int(np.floor((t[0]+180.)/float(res)))
#print >> sys.stderr, xIndx, yIndx, data
outMap[pft,yIndx,xIndx]+=data
outMapN[pft,yIndx,xIndx]+=1.
for p in xrange(int(nPft)):
for y in xrange(int(np.floor((180/float(res))))):
for x in xrange(int(np.floor((360/float(res))))):
#print (x,y,outMapN[y,x]),
if outMapN[p,y,x]>0:
try:
print 0.01*outMap[p,y,x]/outMapN[p,y,x],
except:
print >> sys.stderr, ">>> ",outMap[p,y,x],outMapN[p,y,x]
sys.exit()
else:
print "0",
print
print
if __name__=="__main__":
filename = "MCI2_red_global_sza0_mask.img"
g=globClump( filename )
# map centre
#print "###",g.getPointLatLon( 21600.0, 43200.0 )
#print "###",g.getPointLatLon( 0, 43200.0 )
#print "###",g.getPointLatLon( 0, 0.0 )
#print "###",g.getPointLatLon( 21600.0, 0.0 )
g.getLatLonAsciiMap()
The thing is that the rescaling processes alone with the first map works fine. It also works fine for higher northern latitudes when I merge the two maps, however, it comes with constant strips of zeros in places where it should be the mean like this example of a limited box over Amazon, but all over the world:
