Слияние кода завершено, страница обновится автоматически
/*
* Created by SharpDevelop.
* User: Administrator
* Date: 2016/11/11
* Time: 16:51
*
* To change this template use Tools | Options | Coding | Edit Standard Headers.
*/
using System;
using MapProj;
namespace GeoFly
{
/// <summary>
/// Description of Sampling.
/// </summary>
public static class Sampling
{
private static int Search(double value, double[] array)
{
int len1 = array.Length;
if (value <= array[0])
return 0;
if (value >= array[len1 - 1])
return len1 - 1;
for (int i = 0; i < len1; i++) {
if (value < array[i])
return i - 1;
}
throw new Exception("Invalid result");
}
private static int SearchReverse(double value, double[] array)
{
int len1 = array.Length;
if (value >= array[0])
return 0;
if (value <= array[len1 - 1])
return len1 - 1;
for (int i = 0; i < len1; i++) {
if (value > array[i])
return i - 1;
}
throw new Exception("Invalid result");
}
private static double LinearInter(double AValue, double BValue, double ACoor, double BCoor, double ECoor)
{
if (BCoor - ACoor == 0)
return (AValue + BValue) / 2;
return (AValue * (BCoor - ECoor) + BValue * (ECoor - ACoor)) / (BCoor - ACoor);
}
/// <summary>
/// Interpolate data using bilinear functions
/// newlats,newlons is the destination coordinates
/// oldlats,oldlons is the original coordinates of data
/// example(for the domain of China):
/// ReGrid(data,numpy.arange(0,61),numpy.arange(70,141),lat,lon,ReverseLat=True)
/// </summary>
/// <param name="newlat"></param>
/// <param name="newlon"></param>
/// <returns></returns>
public static double[,] ReGrid_NotEven(double[,] data, double[] latsSource, double[] lonsSource, double[] newlats, double[] newlons, bool ReverseLat)
{
int nlat = newlats.Length;
int nlon = newlons.Length;
double[,] result = new double[nlat, nlon];
int[] lonIndice = new int[nlon];
for (int i = 0; i < nlon; i++) {
double lon = newlons[i];
lonIndice[i] = Search(lon, lonsSource);
}
for (int ilat = 0; ilat < nlat; ilat++) {
double lat = newlats[ilat];
int latIndex = -1;
if (ReverseLat == false) {
latIndex = Search(lat, latsSource);
} else {
latIndex = SearchReverse(lat, latsSource);
}
for (int ilon = 0; ilon < nlon; ilon++) {
//Calculate bilinear interpolation using nearest four grid cells
double lon = newlons[ilon];
//lonIndex=Search(lon,lons1);
int lonIndex = lonIndice[ilon];
double A = data[latIndex, lonIndex];
double B = data[latIndex, lonIndex + 1];
double C = data[latIndex + 1, lonIndex];
double D = data[latIndex + 1, lonIndex + 1];
double lonA = lonsSource[lonIndex];
double lonB = lonsSource[lonIndex + 1];
double E = LinearInter(A, B, lonA, lonB, lon);
double F = LinearInter(C, D, lonA, lonB, lon);
double latA = latsSource[latIndex];
double latB = latsSource[latIndex + 1];
result[ilat, ilon] = LinearInter(E, F, latA, latB, lat);
}
}
return result;
}
public static double[,] ReGridBilinear(double[,] data, double[] latsSource, double[] lonsSource, double[] newlats, double[] newlons)
{
int nlat = newlats.Length;
int nlon = newlons.Length;
double[,] result = new double[nlat, nlon];
int[] lonIndice = new int[nlon];
double dlat = -(latsSource[1] - latsSource[0]);
double dlon = lonsSource[1] - lonsSource[0];
if(dlat<0)
{
Array.Reverse(latsSource);
dlat*=-1;
}
for (int i = 0; i < nlon; i++) {
double lon = newlons[i];
lonIndice[i] = (int)((lon - lonsSource[0]) / dlon);// Search(lon, lonsSource);
}
for (int ilat = 0; ilat < nlat; ilat++) {
double lat = newlats[ilat];
int latIndex = -(int)((lat - latsSource[0]) / dlat);// Search(lat, latsSource);
//latIndex=latsSource.Length-1-latIndex;
for (int ilon = 0; ilon < nlon; ilon++) {
//Calculate bilinear interpolation using nearest four grid cells
double lon = newlons[ilon];
//lonIndex=Search(lon,lons1);
int lonIndex = lonIndice[ilon];
double A = data[latIndex, lonIndex];
double B = A;
if (lonIndex + 1 < lonsSource.Length)
B = data[latIndex, lonIndex + 1];
double C = A;
double D = A;
if (latIndex+1 < latsSource.Length) {
C = data[latIndex+1, lonIndex];
if (lonIndex + 1 < lonsSource.Length)
D = data[latIndex+1, lonIndex + 1];
}
double[] abcd=new double[]{A,B,C,D};
bool flag=false;
for(int ind=0;ind<4;ind++)
{
if(abcd[ind]<-9990)
{
result[nlat- ilat-1,ilon]= Math.Max(A,Math.Max(B,Math.Max(C,D)));
flag=true;
}
}
if(flag)
continue;
double lonA = lonsSource[lonIndex];
double lonB = lonA;
if (lonIndex + 1 < lonsSource.Length)
lonB = lonsSource[lonIndex + 1];
double E = LinearInter(A, B, lonA, lonB, lon);
double F = LinearInter(C, D, lonA, lonB, lon);
double latA = latsSource[latIndex];
double latB = latA;
if (latIndex + 1 < latsSource.Length)
latB = latsSource[latIndex+1];
double vvv= LinearInter(E, F, latA, latB, lat);
result[nlat-1- ilat, ilon] =vvv;
}
}
return result;
}
private static double Interp(double curLat,double curLon,double[,] data,double[] latsSource, double[] lonsSource)
{
double dlat=-(latsSource[1] - latsSource[0]);
double dlon = lonsSource[1] - lonsSource[0];
int latIndex= -(int)((curLat - latsSource[0]) / dlat);
int lonIndex= (int)((curLon - lonsSource[0]) / dlon);
if(dlat<0)
{
dlat*=-1;
}
if(lonIndex+1>=lonsSource.Length)
{
return -9999;
}
double A=data[latIndex ,lonIndex];
double B=data[latIndex,lonIndex+1];
double C=data[latIndex+1,lonIndex];
double D=data[latIndex+1,lonIndex+1];
double lonA=lonsSource[lonIndex];
double lonB=lonsSource[lonIndex+1];
double E=LinearInter(A,B,lonA,lonB,curLon);
double F=LinearInter(C,D,lonA,lonB,curLon);
double latA=latsSource[latIndex];
double latB=latsSource[latIndex+1];
if(A<-9990 || B<-9990 || C<-9990 || D<-9990)
return Math.Max(A,Math.Max(B,Math.Max(C,D)));
return LinearInter(E,F,latA,latB,curLat);
}
public static double[,] ReGridLambert(double[,] data, double[] latsSource, double[] lonsSource, double startLat, double startLon,double cellSize,int rowCount,int colCount
,double refLat,double refLon,double stdLat1,double stdLat2)
{
Array.Reverse(latsSource);
double[,] result = new double[rowCount, colCount];
MapProj.GaussKruger proj=new GaussKruger();
double[] values=proj.ToXY_Lambert(startLon,startLat,refLon,refLat,stdLat1,stdLat2);
double lux=values[0];
double luy=values[1];
double X0=0;
double Y0=0;
for(int row=0;row<rowCount;row++)
{
for(int col=0;col<colCount;col++)
{
X0=lux+cellSize*col;
Y0=luy-cellSize*row;
double[] vv= proj.ToLonLa_Lambert(X0,Y0,refLon,refLat,stdLat1,stdLat2);
double lat= vv[1];
double lon=vv[0];
result[row, col]=Interp(lat,lon,data,latsSource,lonsSource);
}
}
return result;
}
private static double AbsDist(double x0,double y0,double x1,double y1)
{
return Math.Abs(x0-x1)+Math.Abs(y0-y1);
}
public static double[,] ReGridNearest(double[,] data, double[] latsSource, double[] lonsSource, double[] newlats, double[] newlons)
{
int nlat = newlats.Length;
int nlon = newlons.Length;
double[,] result = new double[nlat, nlon];
int[] lonIndice = new int[nlon];
double dlat = -(latsSource[1] - latsSource[0]);
double dlon = lonsSource[1] - lonsSource[0];
if(dlat<0)
{
Array.Reverse(latsSource);
dlat*=-1;
}
for (int i = 0; i < nlon; i++) {
double lon = newlons[i];
lonIndice[i] = (int)((lon - lonsSource[0]) / dlon);// Search(lon, lonsSource);
}
for (int ilat = 0; ilat < nlat; ilat++) {
double lat = newlats[ilat];
int latIndex = -(int)((lat - latsSource[0]) / dlat);// Search(lat, latsSource);
//latIndex=latsSource.Length-1-latIndex;
for (int ilon = 0; ilon < nlon; ilon++) {
//Calculate bilinear interpolation using nearest four grid cells
double lon = newlons[ilon];
//lonIndex=Search(lon,lons1);
int lonIndex = lonIndice[ilon];
double A = data[latIndex, lonIndex];
double ADist=AbsDist(lat,lon,latsSource[latIndex],lonsSource[lonIndex]);
double B = A;
double BDist=ADist;
if (lonIndex + 1 < lonsSource.Length)
{
B = data[latIndex, lonIndex + 1];
BDist=AbsDist(lat,lon,latsSource[latIndex],lonsSource[lonIndex+1]);
}
double C = A;
double CDist=ADist;
double D = A;
double DDist=ADist;
if (latIndex+1 < latsSource.Length) {
C = data[latIndex+1, lonIndex];
CDist= AbsDist(lat,lon,latsSource[latIndex+1],lonsSource[lonIndex]);
if (lonIndex + 1 < lonsSource.Length)
{
D = data[latIndex+1, lonIndex + 1];
DDist=AbsDist(lat,lon,latsSource[latIndex+1],lonsSource[lonIndex+1]);
}
}
double[] aaa=new double[]{A,B,C,D};
double[] bbb=new double[]{ADist,BDist,CDist,DDist};
double min=1e20;
double minValue=0;
for(int i=0;i<4;i++)
{
if(min>bbb[i])
{
min=bbb[i];
minValue=aaa[i];
}
}
result[nlat-1- ilat, ilon] =minValue;
}
}
return result;
}
}
}
Вы можете оставить комментарий после Вход в систему
Неприемлемый контент может быть отображен здесь и не будет показан на странице. Вы можете проверить и изменить его с помощью соответствующей функции редактирования.
Если вы подтверждаете, что содержание не содержит непристойной лексики/перенаправления на рекламу/насилия/вульгарной порнографии/нарушений/пиратства/ложного/незначительного или незаконного контента, связанного с национальными законами и предписаниями, вы можете нажать «Отправить» для подачи апелляции, и мы обработаем ее как можно скорее.
Опубликовать ( 0 )