地理范围内的最大高程与最小高程的获取

liyihongcug

手绘一个范围 需要获取范围内的 最大高程与最小高程，用于湮没分析

--找不到自己独立实现public struct FourP
       {
           public int stx;
           public int xc;
           public int sty;
           public int yc;
           public int zoom;
       } public void PingTu(FourP fp)
      {
          //fp.yc = fp.yc - 1;
          Gdal.AllRegister();
          Driver dr = Gdal.GetDriverByName("GTiff");
          Dataset dst = dr.Create("g:\\baiduZ"+fp.zoom.ToString ()+".tif", fp.xc *256, fp.yc*256, 3, DataType.GDT_Byte, null);
         
          for (int iy = 0; iy < fp.yc; iy++)
          {
              for (int ix = 0; ix < fp.xc; ix++)
              {
                  try
                  {
                      Dataset ds = Gdal.Open("g:\\bdmap\\Z" + fp.zoom.ToString() + "x" + iy.ToString() + "y" + ix.ToString() + ".jpg", Access.GA_ReadOnly);
                      SaveBitmapBuffered(ds, dst, ix, (fp.yc-1)-iy);
                  }
                  catch (Exception ex)
                  {
                      MessageBox.Show(ex.ToString());
                  }
              }          }
          dst.Dispose();
      }  void SaveBitmapBuffered(Dataset ds, Dataset dw, int x, int y)
      {
          Band redBand = ds.GetRasterBand(1);
          Band greenBand = ds.GetRasterBand(2);
          Band blueBand = ds.GetRasterBand(3);
          int width = redBand.XSize;
          int height = redBand.YSize;          DateTime start = DateTime.Now;          byte[] r = new byte[width * height];
          byte[] g = new byte[width * height];
          byte[] b = new byte[width * height];          redBand.ReadRaster(0, 0, width, height, r, width, height, 0, 0);
          greenBand.ReadRaster(0, 0, width, height, g, width, height, 0, 0);
          blueBand.ReadRaster(0, 0, width, height, b, width, height, 0, 0);
          TimeSpan renderTime = DateTime.Now - start;//计算生成一幅新的栅格图像所需要的时间
          listBox1.Items.Add("SaveBitmapBuffered fetch time: " + renderTime.TotalMilliseconds + " ms");
          Band wrb = dw.GetRasterBand(1);
          wrb.WriteRaster(x * width, y * height, width, height, r, width, height, 0, 0);
          Band wgb = dw.GetRasterBand(2);
          wgb.WriteRaster(x * width, y * height, width, height, g, width, height, 0, 0);
          Band wbb = dw.GetRasterBand(3);
          wbb.WriteRaster(x * width, y * height, width, height, b, width, height, 0, 0);          //dst.Dispose();      }
分享：

liyihongcug

import org.gdal.gdal.*;
import org.gdal.ogr.DataSource;
import org.gdal.ogr.Driver;
import org.gdal.ogr.*;
public class VectorToGeojson {
        public static void main(String[] args){
                //注册所有驱动
                ogr.RegisterAll();
                //为了支持中文路径，请添加下面这句代码
                gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
                //为了使属性标志段支持中文，添加下面这句代码
                gdal.SetConfigOption("SHAPE_ENCODING", "");
               
                //使用ogr库打开shape文件
                DataSource ds = ogr.Open("C:\\Users\\Administrator\\Documents\\shpData\\中国地图ArcGIS练习数据\\国界线.shp",0);
                if(ds == null){
                        System.out.println("DataSource is null");
                        return;
                }
                Driver dv = ogr.GetDriverByName("GeoJSON");
                if(dv == null){
                        System.out.println("DataSource is null");
                        return;
                }
                dv.CopyDataSource(ds, "C:\\Users\\Administrator\\Documents\\shpData\\国界线.geojson");
                System.out.println("well done!!");
        }
}

#include "ogrsf_frmts.h"

int main()

{
        OGRRegisterAll();

        OGRDataSource *poDS;

        poDS = OGRSFDriverRegistrar::Open( "G:\\LJF\\point.shp", FALSE );//shape文件存放的路径（point.shp即为自己创建的文件）
        if( poDS == NULL )
        {
                printf( "Open failed.\n%s" );
                exit( 1 );
        }

        OGRLayer  *poLayer;

        poLayer = poDS->GetLayerByName( "point" );

        OGRFeature *poFeature;

        poLayer->ResetReading();
        while( (poFeature = poLayer->GetNextFeature()) != NULL )//获得要素，本实例指的是五个点，所以会循环5次
        {
                OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
                int iField;
                int i=poFDefn->GetFieldCount(); //获得字段的数目，本实例返回5，不包括前两个字段（FID,Shape），这两个字段在arcgis里也不能被修改;
                for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
                {
                        OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );
           //根据字段值得类型，选择对应的输出
                        if( poFieldDefn->GetType() == OFTInteger )
                                printf( "%d,", poFeature->GetFieldAsInteger( iField ) );
                        else if( poFieldDefn->GetType() == OFTReal )
                                printf( "%.3f,", poFeature->GetFieldAsDouble(iField) );
                        else if( poFieldDefn->GetType() == OFTString )
                                printf( "%s,", poFeature->GetFieldAsString(iField) );
                        else
                                printf( "%s,", poFeature->GetFieldAsString(iField) );
                }

                OGRGeometry *poGeometry;

                poGeometry = poFeature->GetGeometryRef();
                if( poGeometry != NULL
                        && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint )
                {
                        OGRPoint *poPoint = (OGRPoint *) poGeometry;

                        printf( "%.3f,%3.f\n", poPoint->getX(), poPoint->getY() );
                }
                else
                {
                        printf( "no point geometry\n" );
                }      
                OGRFeature:estroyFeature( poFeature );
        }

        OGRDataSource::DestroyDataSource( poDS );
        system("pause");
        return 0;
}

liyihongcug

正解 应用GDAL通过传入范围获取tif数据的最大高程值及其坐标
public string GetMultifyElevation(string positions)
        {
            positions = "116.0,40.166667,116.25,40.0";//模拟传入的范围（left，top，right，bottom）
            float dProjX, dProjY , dProjX1, dProjY1;
            string[] arr = positions.Split(',');
            dProjX =float.Parse(arr[0]);
            dProjY = float.Parse(arr[1]);
            dProjX1 = float.Parse(arr[2]);
            dProjY1 = float.Parse(arr[3]);


            string strFilePath = "C:\\webservices\\data\\srtm\\chinaclip.tif";//读取的文件路径
            string testPath = "C:\\webservices\\data\\chinaclip.tif";//要写的文件路径
            string elvate = "";
            Gdal.AllRegister();
            Gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");    //支持中文
            Dataset ds = Gdal.Open(strFilePath, Access.GA_ReadOnly);
            try
            {
                Band Band = ds.GetRasterBand(1);


                //获取图像的尺寸               
                int width = Band.XSize;
                int height = Band.YSize;




                //获取坐标变换系数
                double[] adfGeoTransform = new double[6];
                ds.GetGeoTransform(adfGeoTransform);
                //获取行列号
                double dTemp = adfGeoTransform[1] * adfGeoTransform[5] - adfGeoTransform[2] * adfGeoTransform[4];
                double 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;
                    int dc = Convert.ToInt32(dCol);
                    int dr = Convert.ToInt32(dRow);


                    dCol = (adfGeoTransform[5] * (dProjX1 - adfGeoTransform[0]) -
                           adfGeoTransform[2] * (dProjY1 - adfGeoTransform[3])) / dTemp + 0.5;
                    dRow = (adfGeoTransform[1] * (dProjY1 - adfGeoTransform[3]) -
                        adfGeoTransform[4] * (dProjX1 - adfGeoTransform[0])) / dTemp + 0.5;
                    int dc1 = Convert.ToInt32(dCol);
                    int dr1 = Convert.ToInt32(dRow);
                    int fx = dc - dc1;
                    int fy = dr - dr1;
                    if (fx < 0)
                        fx = -fx;
                    if (fy < 0)
                        fy = -fy;


                    //获取DEM数值到一维数组
                    float[] data = new float[fx * fy];
                    //读取感觉可以去掉
                    CPLErr err = Band.ReadRaster(dc, dr, fx, fy, data, fx, fy, 0, 0);
                    //裁切
                    int[] bandmap = { 1 };
                    DataType DT = DataType.GDT_CFloat32;
                    Dataset dataset = ds.GetDriver().Create(testPath, fx, fy, 1, DT, null);
                    //写入
                    dataset.WriteRaster(0, 0, fx, fy, data, fx, fy, 1, bandmap, 0, 0, 0);


                  


                    Band bd = dataset.GetRasterBand(1);
                    //获取最小最大值
                    double[] lim = new double[2];
                    bd.ComputeRasterMinMax(lim, 1);   //最后的ApproxOK设为1，设为0效果一样
                    float _Min = (float)lim[0];
                    float _Max = (float)lim[1];   


                    bd.ReadRaster(0, 0, fx, fy, data, fx, fy, 0, 0);
                int c =0;
                int x1 = -1, y1 = -1, x2 = -1, y2 = -1;
                //遍历获取行列值
                    for (int i = 0; i < fx; i++)
                    {
                        if (x2 != -1 && y2 != -1 && x1 != -1 && y1 != -1)
                        {
                            goto conmehere;//为了提高效率使用goto语句
                        }


                        for (int j = 0; j < fy; j++)
                        {
                            if (x2 != -1 && y2 != -1 && x1 != -1 && y1 != -1)
                            {
                                goto conmehere;
                            }
                            if (_Min == data[c++])
                            {


                                x1 = i;
                                y1 = j;                              
                            }
                            if (_Max == data[c])
                            {
                                x2 = i;
                                y2 = j;                              
                            }
                        }
                    }
                conmehere:


                    //adfGeoTransform[0]  左上角x坐标   
                    //adfGeoTransform[1]  东西方向分辨率  
                    //adfGeoTransform[2]  旋转角度, 0表示图像 "北方朝上"  
                    //adfGeoTransform[3]  左上角y坐标   
                    //adfGeoTransform[4]  旋转角度, 0表示图像 "北方朝上"  
                    //adfGeoTransform[5]  南北方向分辨率  
                    Band.Dispose();
                    double geox1 = dProjX + x1 * adfGeoTransform[1] + y1 * adfGeoTransform[2];
                    double geoy1 = dProjY + x1 * adfGeoTransform[4] + y1 * adfGeoTransform[5];
                    double geox2 = dProjX + x2 * adfGeoTransform[1] + y2 * adfGeoTransform[2];
                    double geoy2 = dProjY + x2 * adfGeoTransform[4] + y2 * adfGeoTransform[5];


                    elvate = geox1 + "," + geoy1 + "," + _Min + ";" + geox2 + "," + geoy2 + "," + _Max;
                return elvate;
            }
            catch
            {
                return "";
            }
        }问题1、测试程序注册gdal时出现错误

System.MethodAccessException”类型的未经处理的异常出现在 gdal_csharp.dll 中。
其他信息: 安全透明方法“OSGeo.GDAL.Gdal.AllRegister()”尝试通过方法“OSGeo.GDAL.GdalPINVOKE.AllRegister()”调用本机代码失败。方法必须是安全关键的或安全可靠关键的，才能调用本机代码。
解决办法：

修改swig生成的C#封装类代码，强制声明为可被安全透明代码调用即可，以..\swig\csharp\gdal\Gdal.cs类和..\swig\csharp\gdal\Dataset.cs类为例，在其类声明的开头添加下述两行代码:

using System.Security;
[SecuritySafeCritical]