test_gdal_transform.c

#include "gdal.h"
#include "cpl_conv.h" // for CPLMalloc()
#include "gdal_alg.h"
#include <stdio.h>

int main(int argc, char *argv[]) {
    	GDALAllRegister();
	
	//Source Dataset from NetCDF
	GDALDatasetH  hDataset;
	hDataset = GDALOpen("NETCDF:\"/g/data/rr5/satellite/obs/himawari8/FLDK/2016/06/29/0120/20160629012000-P1S-ABOM_BRF_B01-PRJ_GEOS141_2000-HIMAWARI8-AHI.nc\":channel_0001_brf", GA_ReadOnly);

	if(hDataset == NULL) {
		return 1;
	}

	//Destination Dataset as MEM Driver
	float *tile;
	tile = (float *) CPLMalloc(sizeof(float)*256*256);
	char memStr[80];
	sprintf(memStr, "MEM:::DATAPOINTER=%u,PIXELS=%d,LINES=%d,DATATYPE=Float32", &tile[0], 256, 256);
	printf("%s\n", memStr);
	GDALDatasetH  hMDataset;
	hMDataset = GDALOpen(memStr, GA_Update);

	if(hDataset == NULL) {
		return 1;
	}

	//Set Projection and GeoTransform
	GDALSetProjection(hMDataset, "PROJCS[\"WGS 84 / Pseudo-Mercator\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Mercator_1SP\"],PARAMETER[\"central_meridian\",0],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],EXTENSION[\"PROJ4\",\"+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs\"],AUTHORITY[\"EPSG\",\"3857\"]]");
	double destGeot[6] = {12503874.85, 10000, 0, -1408887.35, 0, -10000};
	GDALSetGeoTransform(hMDataset, destGeot);

	//Define Transformer function
	void *trans_cb;
	trans_cb = GDALCreateGenImgProjTransformer(hDataset, GDALGetProjectionRef(hDataset), hMDataset, GDALGetProjectionRef(hMDataset), 0, 0, 0);

	//Call Transform
	int ok = 0;
	double x = 0;
	double y = 0;
	double z = 0;
	GDALGenImgProjTransform(trans_cb, 0, 256*256, &x, &y, &z, &ok);

	CPLFree(tile);
	GDALClose(hDataset);
}