subsetfunc.cpp

提取NOAA PAL(Pathfinder AVHRR Land cover data)中的中国部分子集

#include "stdafx.h"
#include "subsetfunc.h"

#define PAL_SAMPLE		5004
#define PAL_LINE		2168
#define PAL_RESOLUTION	8000.0

#define GOODE_UL_X		-20015500
#define GOODE_UL_Y		8673500

void WriteASCIIHeader(CStdioFile& f, Corner& c, BOOL bNorth);
void WriteEnviHeader(Corner& c, BOOL bNorth);

void SubsetToASCII(CString inputFile, CString outputPath, Corner& nc, Corner& sc, int recordSize, int offset, double gain)
{
	FILE *stream;
	CStdioFile sf1, sf2;
	CFileException fe0, fe1, fe2;

	CString outputFile1 = outputPath + inputFile.Right(inputFile.GetLength()-inputFile.ReverseFind('\\')-1);
	CString outputFile2 = outputFile1+"s.txt";
	outputFile1 += "n.txt";
	if ((stream = fopen(inputFile, "r+b")) != NULL && 
		sf1.Open(outputFile1, CFile::modeCreate | CFile::modeWrite, &fe1) &&
		sf2.Open(outputFile2, CFile::modeCreate | CFile::modeWrite, &fe2))
	{
		WriteASCIIHeader(sf1, nc, TRUE);
		WriteASCIIHeader(sf2, sc, FALSE);

		UINT sample = nc.ULSample;
		UINT line = nc.ULLine;
		long startPos1 = ((line-1)*PAL_SAMPLE+sample)*recordSize;
		int fh = fileno(stream);
		long fileLength = _filelength(fh);
		fseek(stream, startPos1, 0);

		CString temp;
		long pos, count=0;
		while (ftell(stream) <= fileLength)
		{
			if(recordSize == 1)
			{
				int number = fgetc(stream);
				temp.Format("%.3f ", (number-offset)*gain);
			}
			else
			{
				int number1 = fgetc(stream);
				int number2 = fgetc(stream);
				int number3 = (number1<<8)+number2;
				temp.Format("%.3f ", (number3-offset)*gain);
			}
			sf1.WriteString(temp);
			sample++;
			if (sample >= nc.LRSample)
			{
				sample = nc.ULSample;
				line++;
				pos = ((line-1)*PAL_SAMPLE+sample)*recordSize;
				fseek(stream, pos, SEEK_SET);
				temp = "\n";
				sf1.WriteString(temp);
			}
			if (line > nc.LRLine)
				break;
		}

		sample = sc.ULSample;
		line = sc.ULLine;
		startPos1 = ((line-1)*PAL_SAMPLE+sample)*recordSize;
		fh = fileno(stream);
		fileLength = _filelength(fh);
		fseek(stream, startPos1, SEEK_SET);

		while (ftell(stream) <= fileLength)
		{
			if(recordSize == 1)
			{
				int number = fgetc(stream);
				temp.Format("%.3f ", (number-offset)*gain);
			}
			else
			{
				int number1 = fgetc(stream);
				int number2 = fgetc(stream);
				int number3 = (number1<<8)+number2;
				temp.Format("%.3f ", (number3-offset)*gain);
			}
			sf2.WriteString(temp);
			sample++;
			if (sample >= sc.LRSample)
			{
				sample = sc.ULSample;
				line++;
				pos = ((line-1)*PAL_SAMPLE+sample)*recordSize;
				fseek(stream, pos, SEEK_SET);
				temp = "\n";
				sf2.WriteString(temp);
			}
			if (line > sc.LRLine)
				break;
		}

		fclose(stream);
		::DeleteFile(inputFile);
		sf1.Close();
		sf2.Close();
	}
}

void WriteASCIIHeader(CStdioFile& f, Corner& c, BOOL bNorth)
{
	int nCols = c.LRSample-c.ULSample, nRows = c.LRLine-c.ULLine;
	double xllCorner = GOODE_UL_X+c.ULSample*PAL_RESOLUTION,
		yllCorner = GOODE_UL_Y-c.LRLine*PAL_RESOLUTION;
	CString temp;
	temp.Format("ncols         %d\n", nCols);
	f.WriteString(temp);
	temp.Format("nrows         %d\n", nRows);
	f.WriteString(temp);
	temp.Format("xllcorner     %f\n", xllCorner);
	f.WriteString(temp);
	temp.Format("yllcorner     %f\n", yllCorner);
	f.WriteString(temp);
	temp.Format("cellsize      %f\n", PAL_RESOLUTION);
	f.WriteString(temp);
	temp = "NODATA_value  -9999\n";
	f.WriteString(temp);

	//生成投影信息文件
	CString project;
	if(bNorth)
		project = "Projection    MOLLWEIDE\n"
				"Units         METERS\n"
				"Zunits        NO\n"
				"Xshift        -3335846.0\n"
				"Yshift        336410.8\n"
				"Parameters    6370997.0  6370997.0\n"
				"  30  0  0.0 /* Longitude of projection center";
	else
		project = "Projection    SINUSOIDAL\n"
				"Spheroid      SPHERE\n"
				"Units         METERS\n"
				"Zunits        NO\n"
				"Xshift        0.0\n"
				"Yshift        0.0\n"
				"Parameters    \n"
				"6370997.0 /* radius of the sphere of reference\n"
				"  30  0  0.0 /* longitude of central meridian\n"
				"3335846.22854 /* false easting (meters)\n"
				"0.0 /* false northing (meters)";
	CString fileName = f.GetFilePath();
	fileName = fileName.Left(fileName.ReverseFind('.'))+".prj";
	CStdioFile pf;
	pf.Open(fileName, CFile::modeCreate | CFile::modeWrite);
	pf.WriteString(project);
	pf.Close();
}

void SubsetToGrid(CString inputFile, CString outputPath, 
	Corner& nc, Corner& sc, int recordSize, int offset, double gain)
{

}

void SubsetToImg(CString inputFile, CString outputPath, 
	Corner& nc, Corner& sc, int recordSize, int offset, double gain)
{

}

void SubsetToEnvi(CString inputFile, CString outputPath, 
	Corner& nc, Corner& sc, int recordSize, int offset, double gain)
{
	FILE *infs, *outnfs, *outsfs;

	CString outputFile1 = outputPath + inputFile.Right(inputFile.GetLength()-inputFile.ReverseFind('\\')-1);
	CString outputFile2 = outputFile1+"s.img";
	outputFile1 += "n.img";
	if ((infs = fopen(inputFile, "r+b")) != NULL && 
		(outnfs = fopen(outputFile1, "w+b")) != NULL &&
		(outsfs = fopen(outputFile2, "w+b")) != NULL)
	{
		WriteEnviHeader(nc, TRUE);
		WriteEnviHeader(sc, FALSE);

		UINT sample = nc.ULSample;
		UINT line = nc.ULLine;
		long startPos1 = ((line-1)*PAL_SAMPLE+sample)*recordSize;
		int fh = fileno(infs);
		long fileLength = _filelength(fh);
		fseek(infs, startPos1, 0);

		CString temp;
		long pos, count=0;
		while (ftell(infs) <= fileLength)
		{
			int number = fgetc(infs);
			double data = (number-offset)*gain;
			fwrite(&data, sizeof(double), 1, outnfs);
			sample++;
			if (sample >= nc.LRSample)
			{
				sample = nc.ULSample;
				line++;
				pos = ((line-1)*PAL_SAMPLE+sample)*recordSize;
				fseek(infs, pos, SEEK_SET);
			}
			if (line > nc.LRLine)
				break;
		}

		sample = sc.ULSample;
		line = sc.ULLine;
		startPos1 = ((line-1)*PAL_SAMPLE+sample)*recordSize;
		fh = fileno(infs);
		fileLength = _filelength(fh);
		fseek(infs, startPos1, SEEK_SET);

		while (ftell(infs) <= fileLength)
		{
			int number = fgetc(infs);
			double data = (number-offset)*gain;
			fwrite(&data, sizeof(double), 1, outsfs);
			sample++;
			if (sample >= sc.LRSample)
			{
				sample = sc.ULSample;
				line++;
				pos = ((line-1)*PAL_SAMPLE+sample)*recordSize;
				fseek(infs, pos, SEEK_SET);
			}
			if (line > sc.LRLine)
				break;
		}

		fclose(infs);
		::DeleteFile(inputFile);
		fclose(outnfs);
		fclose(outsfs);
	}
}

void WriteEnviHeader(Corner& c, BOOL bNorth)
{

}