We have been using fast_hist and fast_avg to map our MODIS level 1b pixels to a uniform lat/lon grid for plotting. There is one big problem with this approach: the actual ground spacing of a uniform longitude grid changes drastically from the equator to the poles, because longitude lines converge. A better approach is to do the following:
1) pick a map projection
2) define an x,y in grid in meters based on that projection
3) resample the satellite data onto that grid
4) save it as a geotiff file, which maintains all the information about the projection we used
In [1]:
import h5py
from a301lib.geolocate import find_corners
import numpy as np
import pyproj
import pyresample
from pyresample import kd_tree,geometry
from pyresample.plot import area_def2basemap
from matplotlib import pyplot as plt
from a301utils.modismeta_read import parseMeta
from a301utils.a301_readfile import download
data_name='MYD021KM.A2016224.2100.006.2016225153002.h5'
download(data_name)
geom='MYD03.A2016224.2100.006.2016225152335.h5'
download(geom)
index=0
my_name = 'EV_250_Aggr1km_RefSB'
with h5py.File(data_name,'r') as h5_file:
chan1=h5_file['MODIS_SWATH_Type_L1B']['Data Fields'][my_name][index,:,:]
scale=h5_file['MODIS_SWATH_Type_L1B']['Data Fields'][my_name].attrs['reflectance_scales'][...]
offset=h5_file['MODIS_SWATH_Type_L1B']['Data Fields'][my_name].attrs['reflectance_offsets'][...]
chan1_calibrated =(chan1 - offset[index])*scale[index]
with h5py.File(geom) as geo_file:
lon_data=geo_file['MODIS_Swath_Type_GEO']['Geolocation Fields']['Longitude'][...]
lat_data=geo_file['MODIS_Swath_Type_GEO']['Geolocation Fields']['Latitude'][...]
The function is parseMeta and you can treat it as a black box unless you're interested in how regular expressions work in python. It does the same thing that find_corners does, but skips the calculation, since NASA has already computed everything we need more accurately.
In [2]:
corners=parseMeta(data_name)
proj_id = 'laea'
datum = 'WGS84'
lat_0 = '{lat_0:5.2f}'.format_map(corners)
lon_0= '{lon_0:5.2f}'.format_map(corners)
lon_bbox = [corners['min_lon'],corners['max_lon']]
lat_bbox = [corners['min_lat'],corners['max_lat']]
The program that resamples modis data onto a particular projected grid is pyresample. I'll resample the swath using a lambert azimuthal equal area projection
The output will be a 2500 x 2500 array called result which has the channel 1 data resampled onto a grid centered at the lat_0,lon_0 center of the swath. The values will be determined by averaging the nearest neighbors to a particular cell location, using a zone of influence with a radius of 5000 meters, and a kd-tree
The next cell puts the projection into a structure that pyresample understands, and does the resampling
In [3]:
area_dict = dict(datum=datum,lat_0=lat_0,lon_0=lon_0,
proj=proj_id,units='m')
prj=pyproj.Proj(area_dict)
x, y = prj(lon_bbox, lat_bbox)
xsize=2200
ysize=2500
area_id = 'granule'
area_name = 'modis swath 5min granule'
area_extent = (x[0], y[0], x[1], y[1])
area_def = geometry.AreaDefinition(area_id, area_name, proj_id,
area_dict, xsize,ysize, area_extent)
swath_def = geometry.SwathDefinition(lons=lon_data, lats=lat_data)
result = kd_tree.resample_nearest(swath_def, chan1_calibrated.ravel(),
area_def, radius_of_influence=5000, nprocs=2)
print(area_def)
pyresample can take its projection structure and turn it into a basemap instance using area_def2basemap
This works because both pyresample and basemap use the same underlying projection code called pyproj.
In [4]:
plt.close('all')
fig,ax = plt.subplots(1,1, figsize=(8,8))
bmap=area_def2basemap(area_def,ax=ax,resolution='c')
num_meridians=180
num_parallels = 90
vmin=None; vmax=None
col = bmap.imshow(result, origin='upper', vmin=0, vmax=0.4)
label='channel 1'
bmap.drawmeridians(np.arange(-180, 180, num_meridians),labels=[True,False,False,True])
bmap.drawparallels(np.arange(-90, 90, num_parallels),labels=[False,True,True,False])
bmap.drawcoastlines()
fig.colorbar(col, shrink=0.5, pad=0.05).set_label(label)
fig.canvas.draw()
plt.show()
In [5]:
from osgeo import gdal, osr
raster = gdal.GetDriverByName("GTiff")
gformat = gdal.GDT_Float32
channel=result.astype(np.float32)
opacity=0
fill_value=0
g_opts=[]
height,width=result.shape
tiffile='test.tif'
dst_ds = raster.Create(tiffile,width,height,
1,gformat,g_opts)
area=area_def
adfgeotransform = [area.area_extent[0], area.pixel_size_x, 0,
area.area_extent[3], 0, -area.pixel_size_y]
dst_ds.SetGeoTransform(adfgeotransform)
srs = osr.SpatialReference()
srs.ImportFromProj4(area.proj4_string)
srs.SetProjCS(area.proj_id)
srs = srs.ExportToWkt()
dst_ds.SetProjection(srs)
dst_ds.GetRasterBand(1).WriteArray(channel)
del dst_ds
The gdal package has some scripts to dump the details about the tif file and the srs ("spatial reference system")
In [6]:
!gdalinfo test.tif
In [7]:
!gdalsrsinfo test.tif
In [8]:
import gdal
from gdalconst import GA_ReadOnly
data = gdal.Open(tiffile, GA_ReadOnly)
geoTransform = data.GetGeoTransform()
minx = geoTransform[0]
maxy = geoTransform[3]
maxx = minx + geoTransform[1] * data.RasterXSize
miny = maxy + geoTransform[5] * data.RasterYSize
print('\nhere are the projected corners of the x,y raster\n{}'.format([minx, miny, maxx, maxy]))
print('\nhere are the attributes of the data instance:\n{}\n'.format(dir(data)))
data = None
In [ ]:
In [ ]: