/* This example shows how to retrieve lon/lat values from HDF-EOS2 grid data. */

#include <mfhdf.h>
#include <hdf.h>
#include <HdfEosDef.h>
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char **argv)
{
	int32 gridfile1;
	int32 grid1;
	int32 grid1_xdim;
	int32 grid1_ydim;
	float64 grid1_upleft[2];
	float64 grid1_lowright[2];
	int32 grid1_projcode;
	int32 grid1_zone;
	int32 grid1_sphere;
	float64 grid1_params[16];
	int32 grid1_pixreg;
	int32 grid1_origin;
	int32 *grid1_row;
	int32 *grid1_col;
	float64 *grid1_lon;
	float64 *grid1_lat;
	int32 datafield1rank;
	int32 datafield1dimsize[32];
	int32 datafield1type;
	char datafield1dimname[512];
	int16 *datafield1data;
	int32 i, j, k;

	/* Open 'MCD43C4.A2000049.005.2006271205532.hdf' using grid API */
	if ((gridfile1 = GDopen("MCD43C4.A2000049.005.2006271205532.hdf", DFACC_RDONLY)) == -1) {
		fprintf(stderr, "error: cannot open grid 'MCD43C4.A2000049.005.2006271205532.hdf'\n");
		return -1;
	}

	/* Open a grid named 'MCD_CMG_BRDF_0.05Deg' */
	if ((grid1 = GDattach(gridfile1, "MCD_CMG_BRDF_0.05Deg")) == -1) {
		fprintf(stderr, "error: cannot attach to 'MCD_CMG_BRDF_0.05Deg'\n");
		return -1;
	}

	/* Retrieve dimensions and X-Y coordinates of corners from 'MCD_CMG_BRDF_0.05Deg' */
	if ((GDgridinfo(grid1, &grid1_xdim, &grid1_ydim, grid1_upleft, grid1_lowright)) == -1) {
		fprintf(stderr, "error: cannot get grid information from 'MCD_CMG_BRDF_0.05Deg'\n");
		return -1;
	}

	/* Retrieve all GCTP projection information from 'MCD_CMG_BRDF_0.05Deg' */
	if ((GDprojinfo(grid1, &grid1_projcode, &grid1_zone, &grid1_sphere, grid1_params)) == -1) {
		fprintf(stderr, "error: cannot get projection information from 'MCD_CMG_BRDF_0.05Deg'\n");
		return -1;
	}

	/* Retrieve pixel registration information from 'MCD_CMG_BRDF_0.05Deg' */
	if ((GDpixreginfo(grid1, &grid1_pixreg)) == -1) {
		fprintf(stderr, "error: cannot get pixel registration information from 'MCD_CMG_BRDF_0.05Deg'\n");
		return -1;
	}

	/* Retrieve grid pixel origin from 'MCD_CMG_BRDF_0.05Deg' */
	if ((GDorigininfo(grid1, &grid1_origin)) == -1) {
		fprintf(stderr, "error: cannot get pixel origin information from 'MCD_CMG_BRDF_0.05Deg'\n");
		return -1;
	}

	/* Allocate buffer for row */
	if ((grid1_row = malloc(sizeof(int32) * grid1_xdim * grid1_ydim)) == NULL) {
		fprintf(stderr, "error: cannot allocate memory for row\n");
		return -1;
	}
	/* Allocate buffer for column */
	if ((grid1_col = malloc(sizeof(int32) * grid1_xdim * grid1_ydim)) == NULL) {
		fprintf(stderr, "error: cannot allocate memory for column\n");
		return -1;
	}
	/* Fill two arguments, rows and columns */
	for (k = j = 0; j < grid1_ydim; ++j) {
		for (i = 0; i < grid1_xdim; ++i) {
			grid1_row[k] = j;
			grid1_col[k] = i;
			++k;
		}
	}
	/* Allocate buffer for longitude */
	if ((grid1_lon = malloc(sizeof(float64) * grid1_xdim * grid1_ydim)) == NULL) {
		fprintf(stderr, "error: cannot allocate memory for longitude\n");
		return -1;
	}
	/* Allocate buffer for latitude */
	if ((grid1_lat = malloc(sizeof(float64) * grid1_xdim * grid1_ydim)) == NULL) {
		fprintf(stderr, "error: cannot allocate memory for latitude\n");
		return -1;
	}
	/* Retrieve lon/lat values for 'MCD_CMG_BRDF_0.05Deg' */
	if ((GDij2ll(grid1_projcode, grid1_zone, grid1_params, grid1_sphere, grid1_xdim, grid1_ydim, grid1_upleft, grid1_lowright, grid1_xdim * grid1_ydim, grid1_row, grid1_col, grid1_lon, grid1_lat, grid1_pixreg, grid1_origin)) == -1) {
		fprintf(stderr, "error: cannot retrieve lon/lat values for 'MCD_CMG_BRDF_0.05Deg'\n");
		return -1;
	}
	for (i = 0; i < 4; ++i) {
		for (j = 0; j < 5; ++j)
			printf("%5.1lf/%4.1lf, ", grid1_lon[j + grid1_xdim * i], grid1_lat[j + grid1_xdim * i]);
		printf("\n");
	}
	/* Release buffer for row */
	free(grid1_row);
	/* Release buffer for column */
	free(grid1_col);
	/* Release buffer for longitude */
	free(grid1_lon);
	/* Release buffer for latitude */
	free(grid1_lat);

	/* Retrieve information about 'Nadir_Reflectance_Band1' datafield */
	if ((GDfieldinfo(grid1, "Nadir_Reflectance_Band1", &datafield1rank, datafield1dimsize, &datafield1type, datafield1dimname)) == -1) {
		fprintf(stderr, "error: cannot get the field info for 'Nadir_Reflectance_Band1'\n");
		return -1;
	}
	/* Allocate buffer for 'Nadir_Reflectance_Band1' */
	if ((datafield1data = malloc(sizeof(int16) * 3600 * 7200)) == NULL) {
		fprintf(stderr, "error: cannot allocate memory for 'Nadir_Reflectance_Band1'\n");
		return -1;
	}
	/* Read data from 'Nadir_Reflectance_Band1' */
	if ((GDreadfield(grid1, "Nadir_Reflectance_Band1", NULL, NULL, NULL, datafield1data)) == -1) {
		fprintf(stderr, "error: cannot read field 'Nadir_Reflectance_Band1'\n");
		return -1;
	}
	/* Dump data from 'Nadir_Reflectance_Band1' */
	for (i = 0; i < 4; ++i) {
		for (j = 0; j < 5; ++j) {
			printf("%d ", datafield1data[j + 7200 * i]);
		}
		printf("\n");
	}
	/* Release the buffer for 'Nadir_Reflectance_Band1' */
	free(datafield1data);

	/* Close the grid named 'MCD_CMG_BRDF_0.05Deg' */
	if ((GDdetach(grid1)) == -1) {
		fprintf(stderr, "error: cannot detach from 'MCD_CMG_BRDF_0.05Deg'\n");
		return -1;
	}

	/* Close 'MCD43C4.A2000049.005.2006271205532.hdf' */
	if ((GDclose(gridfile1)) == -1) {
		fprintf(stderr, "error: cannot close grid 'MCD43C4.A2000049.005.2006271205532.hdf'\n");
		return -1;
	}

	return 0;
}