#T.Waza, OpenGIS
#26 March 2005
# The current MODIS HDF import (TNTmips v7.0)  georeferences the rasters with only 4 corner control points.
# Also created are  latitude and longitude rasters located in the geolocation folders under the RVC file on import.
# The lat/lon rasters overlay the data rasters and contain a control point for every 4th cell in the data rasters.

# This script creates a new control point georeference for these imported MODIS HDF rasters using the lat/lon rasters.


clear();

class Raster Rin, R_lat, R_lon;
numeric i, j, k, stp;
numeric xsrc, ysrc, zsrc, xdest, ydest, zdest;
numeric lat, lon, numpoints;

GetInputRaster(Rin);
GetInputRaster(R_lat);
GetInputRaster(R_lon);

numeric numlin_in = Rin.$Info.NumLins;
numeric numcol_in = Rin.$Info.NumCols;
# for NASA MODIS HDF, numlin=8120, numcol=5416

numeric numlin_lat = R_lat.$Info.NumLins;
numeric numcol_lat = R_lat.$Info.NumCols;
# for NASA MODIS HDF, numlin=2030, numcol=1354

numeric factor_lin = numlin_in/numlin_lat;
numeric factor_col = numcol_in/numcol_lat;
# for NASA MODIS HDF, factor_lin=4, factor_col=4

# location of lat/lon values in cell. 
# r = 0 ...upper-left corner
#   0.5 ...middle
#     1 ...lower-right corner

numeric r = 0.5;
# for NASA MODIS HDF, Control points are attached at the center of 4 x 4 cells.

#attach georef
array xsrc[3000], ysrc[3000], zsrc[3000];
array xdest[3000], ydest[3000], zdest[3000];

k = 1;
stp = 50;

for i = 1 to numlin_lat step stp				#Line, y-direction
	begin
		for j = 1 to numcol_lat step stp		#Column, x-direction
			begin

				xsrc[k] = factor_col*(j-1+r);
				ysrc[k] = factor_lin*(i-1+r);
				zsrc[k] = 0;

				xdest[k] = R_lon[i,j];
				ydest[k] = R_lat[i,j];
				zdest[k] = 0;

				k = k + 1;
			end
	end

#set number of control points
numpoints = k-1;
print(numpoints);

#now create the georef via dialog
CreateControlPointGeoref(Rin, numpoints, xsrc, ysrc, zsrc, xdest, ydest, zdest);

beep();

CloseRaster(Rin);
CloseRaster(R_lon);
CloseRaster(R_lat);