Gdal gets the specified geo-coordinate cell value, gets the specified geographic range image data

Tag:gdal   gis    data    algorithm   

Gdalimage.h#include "StructDef.h" #include "gdal1.11.2/gdal_priv.h" #include "gdal1.11.2/gdal.h"//#ifdef __ Cplusplus//extern "C" {//#endifstruct Strasterinfo{char filename[255]; Gdaldataset* Pdataset; Xrect<rtsdatatypegeo> rasterrange;int Nbandcount;}; struct Stbitmapinfo{char Filename[255];bool isinters;//whether there are intersecting parts byte* pbuf;int nbandcount;int xstart;int yStart;int xEnd ; int yend;int Xsize;int ysize;stbitmapinfo (): Isinters (False), PBuf (NULL), Nbandcount (0), Xstart (0), Ystart (0), xEnd (0) , Yend (0), xsize (0), ysize (0) {}~stbitmapinfo () {if (PBuf) {free (pBuf);p buf = NULL;}}; static char g_szrasterpath[255] = {0};//Raster image folder full path path static int g_screenwidth = 0;//screen resolution, wide static int g_screenheight = 0; Screen resolution, high//by geographical coordinates get image row number bool Projection2imagerowcol (double *adfgeotransform, double dprojx, double dprojy, int & icol, int &irow);//Get the cell value of the specified point based on geographic location 2015.3.5//cszrasterpath[in]: Image file absolute path Geox, Geoy[in]: location Pbuf[out]: Store the Read cell value, NULL read failure//nbandcount[out]: Number of bands Xsize, Ysize[in]: To get the range of cell values, and Gaumer think 1void getPixvalbygeopos (const char *cszrasterpath, const double GeoX, const double geoy,int *&pbuf, int &nbandcount, int x Size, int ysize);//= 1 = 1//initialization path and resolution void initvariable (const char *szrasterpath, const int screenwidth, const int scree nheight);//Get Image object void Getbuf (const xrect<rtsdatatypegeo> showrange, Stbitmapinfo *&pbitmap, int& nSize) ;//#ifdef __cplusplus//}//#endif

Gdalimage.cpp#include <vector> #include <dirent.h>//#include <io.h> #include <string.h># Include "StructDef.h" #include "DRECT.h" #include "GdalImage.h" #include "gdal1.11.2/gdal_priv.h" #include "gdal1.11.2/ Gdal.h "#include <iostream> #include" mapWndBase.h "//const int dircharmaxlen = 1024; Defines the maximum directory length static vector<strasterinfo> g_vecri;//Initializes all tif,bmp,img information under the specified folder//by geographical coordinates to get the image row number bool Projection2imagerowcol (double *adfgeotransform, double dprojx, double dprojy, int &icol, int &irow) {try{double Dtemp = adfgeotransform[1]*adfgeotransform[5]-adfgeotransform[2]*adfgeotransform[4];d ouble dCol = 0.0, DRow = 0.0;dCol = (adfgeotransform[5]* (dprojx-adfgeotransform[0])-adfgeotransform[2]* (dprojy-adfgeotransform[3]))/dTemp + 0.5; Drow = (adfgeotransform[1]* (dprojy-adfgeotransform[3])-adfgeotransform[4]* (dprojx-adfgeotransform[0])/dTemp + 0.5 ; icol = int (dcol); iRow = int (drow); return true;} catch (...) {return false;}} Gets the cell value of the specified point based on geographic location 2015.3.5//cszRasterpath[in]: Image file absolute path Geox, Geoy[in]: location Pbuf[out]: Store read to cell value, NULL read failure//nbandcount[out]: Number of bands Xsize, Ysize[in]: To get the range of cell values, width and gaumer think 1void getpixvalbygeopos (const char *cszrasterpath, const double GeoX, const double Geoy,int *&pbuf, I NT &nbandcount, int xsize = 1, int ysize = 1) {try{//gdalallregister (); Gdaldataset *pdataset = (gdaldataset*) gdalopen (Cszrasterpath, ga_readonly);d ouble gt[6];//get affine transformation model pdataset-> Getgeotransform (GT);//Get image width and height int rasterxsize = pdataset->getrasterxsize (); int rasterysize = pdataset-> Getrasterysize ();//Get the number of image bands Nbandcount = Pdataset->getrastercount (); int icol, irow;//geolocation to the image row and column position if (! Projection2imagerowcol (GT, Geoy, GeoX, Icol, IRow)) throw 0;//geo-location conversion results if the image range is exceeded if (Icol < 0 | | icol > Rasterxsize | | IRow < 0 | | IRow > Rasterysize) {pBuf = NULL;} Else{pbuf = new Int[nbandcount*xsize*ysize];int *panbandmap= new int[nbandcount];for (int i = 0; i < nbandcount; ++i) p Anbandmap[i] = i+1;//Gets the cell value of each band for the specified row and column position, deposited Pbufpdataset->rasterio (Gf_read, Icol, IRow, Xsize, Ysize, PBuf, Xsize,ysize, Gdt_cint32, Nbandcount, Panbandmap, 4, 0, 0);d elete[] Panbandmap;pa Nbandmap = NULL;}} catch (...) {if (PBUF) delete[] pbuf;pbuf = NULL;} return;} Under Linux use void getfilelist (const char* szpath, vector<string>& vecfile) {dir* DP = Opendir (szpath); if (! dp) RE Turn;if (chdir (szpath) = = 1) return;struct dirent* de;for (de = Readdir (DP); de; de = Readdir (DP)) {if (de-d_ty PE! = Dt_dir)) {string strName (de->d_name); string strformat = Strname.substr (Strname.length () -3, 3); if ("tif" = = Strformat | | "img" ==strformat | | "BMP" ==strformat) Vecfile.push_back (StrName);}} Closedir (DP); return;} Using//void getfilelist under Windows (const char* szpath, vector<string>& vecfile)//{//_finddata_t fdata; Definition file Lookup Structure Object//long Done;//char Tempdir[dircharmaxlen] = {0}; Defines a temporary character array//strcat (TempDir, szpath); Connection string//strcat (TempDir, "\\*.*"); Connection string (search all files)//done = _findfirst (TempDir, &fdata); Start Find file////if (done!-1)//Whether to find success//{int ret = 0;//while (ret! =-1)//define a loop//{//if (fdata.attrib! = _a_subdir)//Determine file properties//{////cout << fdata.name &LT;&L T Endl;//if (strcmp (Fdata.name, "...")!=0 &&//strcmp (Fdata.name, "...")! = 0 &&//strcmp (fdata.name, ".") = 0)//filter//{//char Dir[dircharmaxlen] = {0};//define Character array//strcat (dir, szpath);//Connection Word String//strcat (dir, "\ \"),//strcat (dir, Fdata.name),//vecfile.push_back (String (dir)),//}//}//ret = _findnext (Done, &fdata); Find Next file//}//}//}bool imagerowcol2projection (double *adfgeotransform, int icol, int iRow, double &AMP;DPROJX, double & AMP;DPROJY) {//adfgeotransform[6] array adfgeotransform holds some of the parameters in the affine transformation, respectively meaning see//adfgeotransform[0] upper left corner x coordinate//  ADFGEOTRANSFORM[1] East-West resolution//ADFGEOTRANSFORM[2] rotation angle, 0 means the image "North facing up"//adfgeotransform[3] upper left corner y coordinate//adfgeotransform[4] Rotation angle, 0 means the image "North facing up"//adfgeotransform[5] North-South direction resolution TRY{DPROJX = Adfgeotransform[0] + adfgeotransform[1] * icol + ADFGEOTRANSFO RM[2] * irow;dprojy = adfgeotransform[3] + adfgeotransform[4] * icol + adfgeotransform[5] *Irow;return true;} catch (...) {return false;}} BOOL Getrasterinfo (vector<strasterinfo>& VECRI) {if (!g_szrasterpath) return false;vector<string> Vecfile; Getfilelist (G_szrasterpath, vecfile); Try{//gdalallregister (); for (int i = 0; i < vecfile.size (); ++i) {Gdaldataset *p Dataset = (gdaldataset*) gdalopen (Vecfile[i].c_str (), ga_readonly);d ouble gt[6];//get affine transformation model pdataset-> Getgeotransform (GT);//Get image width and height int rasterxsize = pdataset->getrasterxsize (); int rasterysize = pdataset-> Getrasterysize ();//Get the number of image bands int nbandcount = Pdataset->getrastercount (); Gdalclose (pdataset);pD ataset = Null;double dprojx, dprojy;if (! Imagerowcol2projection (GT, Rasterxsize, Rasterysize, DPROJX, dprojy)) throw 0;strasterinfo RI; Ri.pdataset = NULL; Ri.nbandcount = nbandcount;strcpy (Ri.filename, Vecfile[i].c_str ()); RI.rasterRange.m_left = gt[0]; RI.rasterRange.m_right = DPROJX; RI.rasterRange.m_top = gt[3]; RI.rasterRange.m_bottom = Dprojy;vecri.push_back (RI);}} catch (...) {return false;} Return TRue;} void initvariable (const char *szrasterpath, const int screenwidth, const int screenheight) {if (!szrasterpath | | screenwi DTH <= 0 | | ScreenHeight <= 0) return;strcpy (G_szrasterpath, szrasterpath); g_screenwidth = Screenwidth;g_screenheight = Screenheight;if (!g_vecri.empty ()) g_vecri.clear (); Getrasterinfo (G_VECRI); return;} BOOL Getbitmapinfo (Strasterinfo stri, int* panbandmap, const xrect<rtsdatatypegeo>& showRange,stBitmapInfo& Amp STBITMAP) {int nbandcount;if (7 = = stri.nbandcount) Nbandcount = 3;elsenbandcount = Stri.nbandcount; xrect<rtsdatatypegeo> Rasterrange = stri.rasterrange;try{//If there is a intersecting range, if there is a intersecting part xrect<rtsdatatypegeo> Interrect;if (!interrect.intersectrect (Showrange, Rasterrange)) {strcpy (stbitmap.filename, stri.filename); Stbitmap.isinters = False;if (stri.pdataset) {gdalclose (stri.pdataset); stri.pdataset = NULL;} return true;} if (!stri.pdataset) Stri.pdataset = (gdaldataset*) gdalopen (Stri.filename, ga_readonly); Gdaldataset *pdataset= stri.pdataset;//Get affine transformation Model double gt[6];pD ataset->getgeotransform (GT); int icol_l=0, irow_b=0, icol_r=0, irow_t=0;/ /left, bottom, right, on the cell row number if (! Projection2imagerowcol (GT, Interrect.m_left, Interrect.m_bottom, icol_l, irow_b)) throw 0;if (! Projection2imagerowcol (GT, Interrect.m_right, Interrect.m_top, Icol_r, irow_t)) throw 0;//get row and column width, high int xsize = ABS (ICol_ r-icol_l); int ysize = ABS (irow_b-irow_t);//buffer size int const BUFLEN = nbandcount*g_screenwidth*g_screenheight; BYTE *pbuf = new Byte[buflen];pD Ataset->rasterio (Gf_read, icol_l, irow_t, Xsize, Ysize, PBuf, g_screenwidth,g_screenh  Eight, Gdt_byte, Nbandcount, Panbandmap, Nbandcount, 0, 1);//Finishing Pbufbyte *pbufback;if (1 = = Nbandcount)//single band {//pbufback = New Int[buflen*4];p bufback = (byte*) malloc (sizeof (int) *buflen*4), for (int i = 0; i < Buflen; ++i) {pbufback[i*4] = P Buf[i];p bufback[i*4+1] = Pbuf[i];p bufback[i*4+2] = Pbuf[i];p bufback[i*4+3] = 0;}} else if (3==nbandcount | | 7==nbandcount)//3 band {//pbufback = new INT[BUFLEN+BUFLEN/3];p Bufback = (byte*) malloc (sizeof (BYTE) * (BUFLEN+BUFLEN/3)); for (int i = 0; i < BUFLEN/3; ++i) {pbufback[i*4] = pbuf[i*3];p bufback[i* 4+1] = pbuf[i*3+1];p bufback[i*4+2] = pbuf[i*3+2];p bufback[i*4+3] = 0;}} else if (4 = = Nbandcount)//4 Band {//pbufback = new Int[buflen];p bufback = (byte*) malloc (buflen*sizeof (int)); for (int i = 0; i < Buflen; ++i) {if (3 = = (i%4)) pbufback[i] = 0;}} Delete[] PBuf; Xrect<rtsdatatypedevice> Bmpdevice; Geotodevice (Interrect,bmpdevice); strcpy (Stbitmap.filename, stri.filename); stbitmap.isinters = true; Stbitmap.nbandcount = Nbandcount;stbitmap.pbuf = pbufback;stbitmap.xstart=bmpdevice.m_left;stbitmap.ystart= Bmpdevice.m_top;stbitmap.xend=bmpdevice.m_right;stbitmap.yend=bmpdevice.m_bottom;stbitmap.xsize = G_screenWidth; Stbitmap.ysize = G_screenheight;} catch (...) {return false;} return true;} void Deletebmparray (Stbitmapinfo *&pbitmap) {delete[] pbitmap;} void Getbuf (const xrect<rtsdatatypegeo> showrange, Stbitmapinfo *&pbitmap, int& nBitmap) {nBitmap = 0;int nbandcount;//band int *panbandmap;int NVECRI = g_vecri.size (); Pbitmap = new Stbitmapinfo[nvecri];for (int i = 0; i < NVECRI; ++i) {nbandcount = G_VECRI[I].NBANDCOUNT;//1 band or 3 band if ( 1==nbandcount | | 3==nbandcount) {panbandmap = new int[nbandcount];for (int j = 0; j < Nbandcount; ++j) panbandmap[j] = j+1;//nbandcount- J;} If 7-band else if (7==nbandcount) {nbandcount = 3;panbandmap = new Int[nbandcount];p anbandmap[0] = 3;panbandmap[1] = 4;panba NDMAP[2] = 2;} if (Getbitmapinfo (G_vecri[i], Panbandmap, Showrange, Pbitmap[nbitmap])) {//Open only the first 4 if ((Pbitmap[nbitmap].isinters) can be displayed && (++nbitmap==4)) {//scan the remaining rasterinfo, turn off open pdatasetfor (int j = i+1; j < Nvecri; ++j) {if (G_vecri[j].pdataset) {gdalclose (g_vecri[j].pdataset); g_vecri[j].pdataset = NULL;}} Break;}} Delete[] Panbandmap;} return;}

Gdal gets the specified geo-coordinate cell value, gets the specified geographic range image data