Geoloc correction for HDF data using gdal

In the previous blog, I explained how to use the gdalwarp tool provided by gdal for correction. In essence, it should be the same as the GLT correction of ENVI. I encapsulated the gdalwarp code and wrote a class for geoloc processing. I hope it will be useful to you.

First, the header file. The comments of the function are very detailed. The following source file is from gdalwarp. if you are interested, see gdalwarp. CPP, the following code just deletes the code that is not used in this file, and some parameters are directly written to death. You can directly look at the original code.

/*************************************** ************************************* Time: * Project: remote sensing processing algorithm * purpose: Geographic lookup table correction * Author: Li minlu * copyright (c) 2013, liminlu0314@gmail.com * describe: Geographic lookup table correction, ********************************** **************************************** * **/# ifndef imagegeolocwarp_h # define imagegeolocwarp_h # include "imgcore. H "/*** \ file imagegeolocwarp. H * @ B Rief **, which is mainly used for correction of HDF data. For example, the marine satellite data and the modemdata. * Generally, this type of data contains two 32-bit floating point data sets of latitude and longpolling. * It is used to store the longitude and latitude coordinates of the data, use these two methods in the correction process */class imgalg_api cimagegeolocwarp {public:/*** @ brief constructor, geolocation table correction * @ Param pszsrcfile raw data path * @ Param pszdstfile result data path * @ Param pszformat result data format * @ Param pprocess progress bar pointer */cimagegeolocwarp (const char * pszsrcfile, const char * pszdstfile, const char * pszformat = "HFA", cprocessbase * pprocess = NULL);/*** @ brief destructor */~ Cimagegeolocwarp (); /*** @ brief get output parameter information * @ Param iheight output image height * @ Param iwidth output image width * @ Param padfgeotransform output geotransform six parameters * @ Param padfextent output image four to the range * @ return: whether the calculation is successful, and Error Code */INT getsuggestedwarpoutput (Int & iheight, Int & iwidth, double * padfgeotransform = NULL, double * padfextent = NULL ); /*** @ brief executes the calibration algorithm * @ Param dresx horizontal resolution of the output image. The default value is 0, indicating that the image size is automatically calculated. * @ Param dresy vertical resolution of the output image, the default value is 0, indicating that the output metadata size is automatically calculated * @ Param RESA Mpling re-sampling method. The default value is bilinear interpolation * @ return, and the error code */INT dogeolocwarp (double dresx = 0.0, double dresy = 0.0, resamplemethod resampling = rm_bilinear); Private: // enter the const char * m_pszsrcfile parameter ;/*! Raw data path */const char * m_pszdstfile ;/*! Result data path */const char * m_pszformat ;/*! Result data format */cproccorner * m_pprocess ;/*! Progress bar pointer */resamplemethod m_resampling ;/*! Sampling Method */string m_strwkt ;/*! Output file projection string, which can only be the longitude and latitude of WGS84 */}; # endif // imagegeolocwarp_h

The following is the source file, which is used to streamline gdalwarp. cpp and remove code unrelated to geoloc.

/*************************************** ************************************* Time: * Project: remote sensing processing algorithm * purpose: Geographic lookup table correction * Author: Li minlu * copyright (c) 2013, liminlu0314@gmail.com * describe: Geographic lookup table correction, ********************************** **************************************** * **/# include "imagegeolocwarp. H "// the header file # include" imagecommon. H "// This header file is the previous blog that uses gdal to manage data (the function used is to delete temporary files )# Include "gdal_priv.h" # include "gdalwarper. H "# include" character "character: cimagegeolocwarp (const char * pszsrcfile, const char * pszdstfile, const char * pszformat, cprocessbase * pprocess) {m_pszsrcfile = pszsrcfile ;/*! Raw data path */m_pszdstfile = pszdstfile ;/*! Result data path */m_pszformat = pszformat ;/*! Result data format */m_pprocess = pprocess ;/*! Progress bar pointer */m_resampling = rm_bilinear ;/*! The default value is bilinear */m_strwkt = srs_wkt_wgs84. // The projection value can only be WGS84} cimagegeolocwarp ::~ Cimagegeolocwarp () {}int cimagegeolocwarp: getsuggestedwarpoutput (Int & iheight, Int & iwidth, double * padfgeotransform, double * padfextent) {If (m_pszsrcfile = NULL | m_pszdstfile = NULL) return re_paramerror; gdalallregister (); gdaldataseth hsrcds = gdalopen (m_pszsrcfile, ga_readonly ); // open the file if (hsrcds = NULL) {If (m_pprocess! = NULL) m_pprocess-> setmessage ("opening the image to be corrected failed! "); Return re_filenotexist;} If (m_strwkt.empty () // The target projection does not exist {If (m_pprocess! = NULL) m_pprocess-> setmessage ("the target projection does not exist! "); Gdalclose (hsrcds); Return re_paramerror;} // create the coordinate conversion parameter char ** papszwarpoptions = NULL; papszwarpoptions = conditions (papszwarpoptions," method "," geoloc_array "); papszwarpoptions = outputs (outputs, "dst_srs", outputs (); void * htransformarg = outputs (hsrcds, null, papszwarpoptions); If (htransformarg = NULL) {If (m_pprocess! = NULL) m_pprocess-> setmessage ("coordinate conversion parameter error! "); Csldestroy (papszwarpoptions); gdalclose (hsrcds); Return re_paramerror;} gdaltransformerinfo * psinfo = (gdaltransformerinfo *) htransformarg; // obtain the target file size and conversion parameter information. Double adfdstgeotransform [6] = {0}; double adfextent [4]; int npixels = 0, nlines = 0; cplerr eerr = convert (hsrcds, psinfo-> pfntransform, htransformarg, adfdstgeotransform, & npixels, & nlines, adfextent, 0); iwidth = pixnels; I Height = nlines; If (padfgeotransform! = NULL) memcpy (padfgeotransform, adfdstgeotransform, sizeof (double) * 6); If (padfextent! = NULL) memcpy (padfextent, adfextent, sizeof (double) * 4); transform (htransformarg); csldestroy (papszwarpoptions); gdalclose (hsrcds); If (eerr! = Ce_none) return re_paramerror; elsereturn re_success;} int cimagegeolocwarp: dogeolocwarp (double dresx, double dresy, resamplemethod resampling) {m_resampling = resampling; if (m_pszsrcfile = NULL | m_pszdstfile = NULL) return re_paramerror; If (fileisused (m_pszdstfile, m_pprocess) return re_fileisuesed; If (m_pprocess! = NULL) {m_pprocess-> resetprocess (); m_pprocess-> setmessage ("image geographic correction is being executed... ");} gdalallregister (); gdaldataseth hsrcds = gdalopen (m_pszsrcfile, ga_readonly); // open the file if (hsrcds = NULL) {If (m_pprocess! = NULL) m_pprocess-> setmessage ("opening the image to be corrected failed! "); Return re_filenotexist;} If (m_strwkt.empty () // The target projection does not exist {If (m_pprocess! = NULL) m_pprocess-> setmessage ("the target projection does not exist! "); Gdalclose (hsrcds); Return re_paramerror;} // create the coordinate conversion parameter char ** papszwarpoptions = NULL; papszwarpoptions = conditions (papszwarpoptions," method "," geoloc_array "); papszwarpoptions = outputs (outputs, "dst_srs", outputs (); void * htransformarg = outputs (hsrcds, null, papszwarpoptions); If (htransformarg = NULL) {If (m_pprocess! = NULL) m_pprocess-> setmessage ("coordinate conversion parameter error! "); Csldestroy (papszwarpoptions); gdalclose (hsrcds); Return re_paramerror;} gdaltransformerinfo * psinfo = (gdaltransformerinfo *) htransformarg; // obtain the target file size and conversion parameter information. Double adfdstgeotransform [6] = {0}; double adfextent [4]; int npixels = 0, nlines = 0; cplerr eerr = gdalsuggestedwarpoutput2 (hsrcds, psinfo-> pfntransform, htransformarg, adfdstgeotransform, & npixels, & nlines, adfext, 0); If (eerr! = Ce_none) {If (m_pprocess! = NULL) m_pprocess-> setmessage ("An error occurred while obtaining the target file size and conversion parameter information! "); Gdaldestroygenimgprojtransformer (htransformarg); csldestroy (papszwarpoptions); gdalclose (hsrcds); Return re_paramerror;} // If you specify the output image resolution if (dresx! = 0.0 | dresy! = 0.0) {// if only one is specified, use the automatically calculated result if (dresx = 0.0) dresx = adfdstgeotransform [1]; If (dresy = 0.0) dresy = adfdstgeotransform [5]; // ensure the resolution is correct if (dresx <0.0) dresx =-dresx; If (dresy> 0.0) dresy =-dresy; // calculate the output image range. Double Minx = adfdstgeotransform [0]; double Maxx = adfdstgeotransform [0] + adfdstgeotransform [1] * npixels; double Maxy = adfdstgeotransform [3]; double miny = adfdstgeotransform [3] + Adfdstgeotransform [5] * nlines; // calculate the output size and range of the image according to the resolution specified by the user. npixels = (INT) (Maxx-Minx)/dresx) + 0.5); nlines = (INT) (miny-Maxy)/dresy) + 0.5); adfdstgeotransform [0] = Minx; adfdstgeotransform [3] = Maxy; adfdstgeotransform [1] = dresx; adfdstgeotransform [5] = dresy;} // obtain the band number int nbandcount = gdalgetrastercount (hsrcds ); // create an output file with the same data type as the source file gdaldatatype EDT = gdalgetrasterdatatype (gdalg Etrasterband (hsrcds, 1); // create the output result file gdaldriverh hdriver = gdalgetdriverbyname (m_pszformat); If (hdriver = NULL) {If (m_pprocess! = NULL) m_pprocess-> setmessage ("this type of data is not supported! "); Round (htransformarg); csldestroy (strong); gdalclose (hsrcds); Return re_filenotsupport;} gdaldataseth hdstds = gdalcreate (hdriver, train, pixnels, nlines, nbandcount, EDT, null); If (hdstds = NULL) {If (m_pprocess! = NULL) m_pprocess-> setmessage ("An error occurred while creating the output file! "); Transform (htransformarg); csldestroy (papszwarpoptions); gdalclose (hsrcds); Return re_createfailed;} // set the information such as projection gdalsetprojection (hdstds, projection ()); gdalsetgeotransform (hdstds, adfdstgeotransform); If (htransformarg! = NULL) gdalsetgenimgprojtransformerdstgeotransform (htransformarg, adfdstgeotransform); // obtain the color table gdalcolortableh HCT one by one; for (INT I = 1; I <= nbandcount; I ++) // band starts from 1 to {HCT = gdalgetrastercolortable (gdalgetrasterband (hsrcds, I); If (Hct! = NULL) gdalsetrastercolortable (gdalgetrasterband (hdstds, I), HCT);} // creates the conversion option csldestroy (papszwarpoptions); // first releases the original papszwarpoptions = NULL; papszwarpoptions = conditions (papszwarpoptions, "init_dest", "0"); optional pfntransformer = NULL; pfntransformer = conditions; // Coordinate Conversion Function gdalwarpoptions * pswarpoptions = conditions (); pswarpoptions-> papszwarpoptions = csldup Licate (random); pswarpoptions-> hsrcds = hsrcds; pswarpoptions-> hdstds = hdstds; pswarpoptions-> nbandcount = nbandcount; // apply for memory pswarpoptions-> random = (int *) cplmalloc (sizeof (INT) * pswarpoptions-> nbandcount); pswarpoptions-> pandstbands = (int *) cplmalloc (sizeof (INT) * pswarpoptions-> nbandcount ); for (Int J = 0; j <nbandcount; j ++) {pswarpoptions-> pansrcbands [J] = J + 1; pswarpoptions-> P Andstbands [J] = J + 1;} pswarpoptions-> progress = (Progress) resampling; pswarpoptions-> pfnprogress = progress; // progress bar callback function; pswarpoptions-> pprogressarg = m_pprocess; // progress bar pointer pswarpoptions-> conditions = conditions; pswarpoptions-> pfntransformer = pfntransformer; conditions-> ptransformerarg = condition; cplerr err = ce_none; condition ooperation; If (ooperation. initi Alize (pswarpoptions) = ce_none) Err = ooperation. values (0, 0, npixels, nlines); else {csldestroy (papszwarpoptions); values (htransformarg); gdalclose (hdstds); gdalclose (hsrcds); rasterdelete (m_pszdstfile ); // return re_failed;} csldestroy (papszwarpoptions); gdaldestroygenimgprojtransformer (htransformarg); gdalclose (hdstds); gdalclose (hsrcds); If (M _ Pprocess &&! M_pprocess-> m_biscontinue) {rasterdelete (m_pszdstfile); // return re_usercancel;} If (Err! = Ce_none) return re_failed; return re_success ;}

The call code is as follows:

# Include <stdio. h> # include <stdlib. h> # include <boost/progress. HPP> // boost timer function int main () {int irev = re_success; cconsoleprocess * ppro = new cconsoleprocess (); progress_timer * ptime = new progress_timer (); // start timing // process a word data in HDF data. The following are the path of the sub-dataset and the path of the output file const char * pszsrc = "hdf4_sds: Unknown: \ "D: \ data \ HDF \ h1bclr101106020418636. l2c. HDF \ ": 12"; const char * pszdst = "d :\\ data \ HDF \ h1bclr101106020418636_1 2. TIF "; cimagegeolocwarp warp (pszsrc, pszdst," HFA ", ppro); irev = warp. dogeolocwarp (0.0, 0.0, rm_nearestneighbour); If (irev = re_success) printf ("success! \ N "); elseprintf (" failed! \ N "); Delete ppro; Delete ptime;: System (" pause "); Return 0 ;}

Finally, the return value in the code, which starts with "RE _", is searched in my previous blog, which is a constant used to mark the returned information. The progress bar class is also described in my previous blog. The third-party libraries used in the above Code are boost and gdal. There is only one gdal In the algorithm, and the boost library is only used for timing in the main function. Other blogs have explained this.

Don't expect that my code can be compiled directly in the past. In that case, I don't think there is much progress for you. If a person understands my code, the functions he doesn't know can probably guess what it means. It's easy to write one by yourself. To learn other people's code, do not expect it to be used directly. The key is to look at other people's ideas, and then use their own projects (projects) for transformation to achieve their own goals, this is a learning method that is advantageous to you. There is a lot of nonsense. I don't think there are any of these things behind it.