KJ135
/
ISCE_INSAR
mirror of http://172.16.0.12:5000/WBJY/ISCE_INSAR.git
4
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Bump autoRIFT 1.0.6 -> 1.0.7 (#191) 

Co-authored-by: Ryan Burns <rtburns-jpl@users.noreply.github.com>
LT1AB
Ryan Burns 2020-09-14 12:52:11 -07:00 committed by GitHub
parent 337dc9c158
commit ba2cb412ea
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 669 additions and 479 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
contrib/geo_autoRIFT/README.md
View File

@ -3,9 +3,9 @@
[![Language](https://img.shields.io/badge/python-3.6%2B-blue.svg)](https://www.python.org/)
[![Latest version](https://img.shields.io/badge/latest%20version-v1.0.6-yellowgreen.svg)](https://github.com/leiyangleon/autoRIFT/releases)
[![Latest version](https://img.shields.io/badge/latest%20version-v1.0.7-yellowgreen.svg)](https://github.com/leiyangleon/autoRIFT/releases)
[![License](https://img.shields.io/badge/License-Apache%202.0-blue.svg)](https://github.com/leiyangleon/autoRIFT/blob/master/LICENSE)
[![Citation](https://img.shields.io/badge/DOI-10.5281/zenodo.3756192-blue)](https://doi.org/10.5281/zenodo.3756192)
[![Citation](https://img.shields.io/badge/DOI-10.5281/zenodo.4025445-blue)](https://doi.org/10.5281/zenodo.4025445)
@ -22,7 +22,7 @@ Copyright (C) 2019 California Institute of Technology. Government Sponsorship A
Link: https://github.com/leiyangleon/autoRIFT
Citation: https://doi.org/10.5281/zenodo.3756192
Citation: https://doi.org/10.5281/zenodo.4025445
## 1. Authors

1
contrib/geo_autoRIFT/SConscript
View File

@ -22,6 +22,7 @@ import os
Import('envcontrib')
package = 'geo_autoRIFT'
envgeoAutorift = envcontrib.Clone()
envgeoAutorift.MergeFlags('-std=c++11')
envgeoAutorift['PACKAGE'] = envcontrib['PACKAGE'] + '/' + package
install = envcontrib['PRJ_SCONS_INSTALL'] + '/' + envgeoAutorift['PACKAGE']
listFiles = ['__init__.py']

1
contrib/geo_autoRIFT/autoRIFT/bindings/autoriftcoremodule.cpp
View File

@ -41,6 +41,7 @@
#include "iostream"
#include "numpy/arrayobject.h"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/imgproc/types_c.h"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/core/core.hpp"

130
contrib/geo_autoRIFT/geogrid/GeogridOptical.py
View File

@ -155,10 +155,10 @@ class GeogridOptical():
# For now print inputs that were obtained
print("Optical Image parameters: ")
print("\nOptical Image parameters: ")
print("X-direction coordinate: " + str(self.startingX) + " " + str(self.XSize))
print("Y-direction coordinate: " + str(self.startingY) + " " + str(self.YSize))
print("Dimensions: " + str(self.numberOfSamples) + " " + str(self.numberOfLines))
print("Dimensions: " + str(self.numberOfSamples) + " " + str(self.numberOfLines) + "\n")
print("Map inputs: ")
print("EPSG: " + str(self.epsgDem))
@ -167,25 +167,45 @@ class GeogridOptical():
print("XLimits: " + str(self._xlim[0]) + " " + str(self._xlim[1]))
print("YLimits: " + str(self._ylim[0]) + " " + str(self._ylim[1]))
print("Extent in km: " + str((self._xlim[1]-self._xlim[0])/1000.0) + " " + str((self._ylim[1]-self._ylim[0])/1000.0))
if (self.demname != ""):
print("DEM: " + str(self.demname))
print("Velocities: " + str(self.vxname) + " " + str(self.vyname))
print("Search Range: " + str(self.srxname) + " " + str(self.sryname))
print("Chip Size Min: " + str(self.csminxname) + " " + str(self.csminyname))
print("Chip Size Max: " + str(self.csmaxxname) + " " + str(self.csmaxyname))
print("Stable Surface Mask: " + str(self.ssmname))
if (self.dhdxname != ""):
print("Slopes: " + str(self.dhdxname) + " " + str(self.dhdyname))
print("Output Nodata Value: " + str(self.nodata_out))
if (self.vxname != ""):
print("Velocities: " + str(self.vxname) + " " + str(self.vyname))
if (self.srxname != ""):
print("Search Range: " + str(self.srxname) + " " + str(self.sryname))
if (self.csminxname != ""):
print("Chip Size Min: " + str(self.csminxname) + " " + str(self.csminyname))
if (self.csmaxxname != ""):
print("Chip Size Max: " + str(self.csmaxxname) + " " + str(self.csmaxyname))
if (self.ssmname != ""):
print("Stable Surface Mask: " + str(self.ssmname))
print("\nOutputs: ")
print("Outputs: ")
print("Window locations: " + str(self.winlocname))
if (self.dhdxname != ""):
if (self.vxname != ""):
print("Window offsets: " + str(self.winoffname))
print("Window search range: " + str(self.winsrname))
print("Window chip size min: " + str(self.wincsminname))
print("Window chip size max: " + str(self.wincsmaxname))
print("Window stable surface mask: " + str(self.winssmname))
print("Window rdr_off2vel_x vector: " + str(self.winro2vxname))
print("Window rdr_off2vel_y vector: " + str(self.winro2vyname))
if (self.srxname != ""):
print("Window search range: " + str(self.winsrname))
if (self.csminxname != ""):
print("Window chip size min: " + str(self.wincsminname))
if (self.csmaxxname != ""):
print("Window chip size max: " + str(self.wincsmaxname))
if (self.ssmname != ""):
print("Window stable surface mask: " + str(self.winssmname))
print("Output Nodata Value: " + str(self.nodata_out) + "\n")
print("Starting processing .... ")
@ -277,6 +297,8 @@ class GeogridOptical():
print("Ylimits : " + str(geoTrans[3] + (lOff + lCount) * geoTrans[5]) + " " + str(geoTrans[3] + lOff * geoTrans[5]))
print("Origin index (in DEM) of geogrid: " + str(pOff) + " " + str(lOff))
print("Dimensions of geogrid: " + str(pCount) + " x " + str(lCount))
projDem = osr.SpatialReference()
@ -296,12 +318,14 @@ class GeogridOptical():
if (self.vxname != ""):
nodata = vxDS.GetRasterBand(1).GetNoDataValue()
else:
nodata = 0
nodata_out = self.nodata_out
pszFormat = "GTiff"
adfGeoTransform = ( self._xlim[0], (self._xlim[1]-self._xlim[0])/pCount, 0, self._ylim[1], 0, (self._ylim[0]-self._ylim[1])/lCount )
adfGeoTransform = ( geoTrans[0] + pOff * geoTrans[1], geoTrans[1], 0, geoTrans[3] + lOff * geoTrans[5], 0, geoTrans[5] )
oSRS = osr.SpatialReference()
pszSRS_WKT = projDem.ExportToWkt()
@ -323,7 +347,7 @@ class GeogridOptical():
if ((self.dhdxname != "")&(self.vxname != "")):
poDriverOff = gdal.GetDriverByName(pszFormat)
if( poDriverOff is None ):
raise Exception('Cannot create gdal driver for output')
@ -339,7 +363,7 @@ class GeogridOptical():
poBand2Off.SetNoDataValue(nodata_out)
if ((self.dhdxname != "")&(self.srxname != "")):
poDriverSch = gdal.GetDriverByName(pszFormat)
if( poDriverSch is None ):
raise Exception('Cannot create gdal driver for output')
@ -354,7 +378,7 @@ class GeogridOptical():
poBand1Sch.SetNoDataValue(nodata_out)
poBand2Sch.SetNoDataValue(nodata_out)
if (self.csminxname != ""):
poDriverMin = gdal.GetDriverByName(pszFormat)
if( poDriverMin is None ):
raise Exception('Cannot create gdal driver for output')
@ -369,7 +393,7 @@ class GeogridOptical():
poBand1Min.SetNoDataValue(nodata_out)
poBand2Min.SetNoDataValue(nodata_out)
if (self.csmaxxname != ""):
poDriverMax = gdal.GetDriverByName(pszFormat)
if( poDriverMax is None ):
raise Exception('Cannot create gdal driver for output')
@ -385,7 +409,7 @@ class GeogridOptical():
poBand2Max.SetNoDataValue(nodata_out)
if (self.ssmname != ""):
poDriverMsk = gdal.GetDriverByName(pszFormat)
if( poDriverMsk is None ):
raise Exception('Cannot create gdal driver for output')
@ -401,7 +425,7 @@ class GeogridOptical():
if (self.dhdxname != ""):
poDriverRO2VX = gdal.GetDriverByName(pszFormat)
if( poDriverRO2VX is None ):
raise Exception('Cannot create gdal driver for output')
@ -417,8 +441,6 @@ class GeogridOptical():
poBand2RO2VX.SetNoDataValue(nodata_out)
poDriverRO2VY = gdal.GetDriverByName(pszFormat)
if( poDriverRO2VY is None ):
raise Exception('Cannot create gdal driver for output')
@ -512,6 +534,8 @@ class GeogridOptical():
slp = np.array([sxLine[jj], syLine[jj], -1.0])
if (self.vxname != ""):
vel = np.array([vxLine[jj], vyLine[jj], 0.0])
else:
vel = np.array([0., 0., 0.])
if (self.srxname != ""):
schrng1 = np.array([srxLine[jj], sryLine[jj], 0.0])
schrng2 = np.array([-srxLine[jj], sryLine[jj], 0.0])
@ -576,22 +600,33 @@ class GeogridOptical():
# else:
# raster11[jj] = 0.
# raster22[jj] = 0.
if (self.vxname == ""):
# pdb.set_trace()
raster11[jj] = 0.
raster22[jj] = 0.
elif (vel[0] == nodata):
if (self.dhdxname != ""):
if (self.vxname != ""):
if (vel[0] == nodata):
raster11[jj] = 0.
raster22[jj] = 0.
else:
raster11[jj] = np.round(np.dot(vel,xunit)*self.repeatTime/self.XSize/365.0/24.0/3600.0*1)
raster22[jj] = np.round(np.dot(vel,yunit)*self.repeatTime/self.YSize/365.0/24.0/3600.0*1)
cross = np.cross(xunit,yunit)
cross = cross / np.linalg.norm(cross)
cross_check = np.abs(np.arccos(np.dot(normal,cross))/np.pi*180.0-90.0)
if (self.srxname == ""):
sr_raster11[jj] = nodata_out
sr_raster22[jj] = nodata_out
elif (vel[0] == nodata):
if (cross_check > 1.0):
raster1a[jj] = normal[2]/(self.repeatTime/self.XSize/365.0/24.0/3600.0)*(normal[2]*yunit[1]-normal[1]*yunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2]));
raster1b[jj] = -normal[2]/(self.repeatTime/self.YSize/365.0/24.0/3600.0)*(normal[2]*xunit[1]-normal[1]*xunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2]));
raster2a[jj] = -normal[2]/(self.repeatTime/self.XSize/365.0/24.0/3600.0)*(normal[2]*yunit[0]-normal[0]*yunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2]));
raster2b[jj] = normal[2]/(self.repeatTime/self.YSize/365.0/24.0/3600.0)*(normal[2]*xunit[0]-normal[0]*xunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2]));
else:
raster1a[jj] = nodata_out
raster1b[jj] = nodata_out
raster2a[jj] = nodata_out
raster2b[jj] = nodata_out
if (self.srxname != ""):
if ((self.vxname != "")&(vel[0] == nodata)):
sr_raster11[jj] = 0
sr_raster22[jj] = 0
else:
@ -609,44 +644,41 @@ class GeogridOptical():
if (self.csminxname != ""):
csmin_raster11[jj] = csminxLine[jj] / self.chipSizeX0 * ChipSizeX0_PIX_X
csmin_raster22[jj] = csminyLine[jj] / self.chipSizeX0 * ChipSizeX0_PIX_Y
else:
csmin_raster11[jj] = nodata_out
csmin_raster22[jj] = nodata_out
if (self.csmaxxname != ""):
csmax_raster11[jj] = csmaxxLine[jj] / self.chipSizeX0 * ChipSizeX0_PIX_X
csmax_raster22[jj] = csmaxyLine[jj] / self.chipSizeX0 * ChipSizeX0_PIX_Y
else:
csmax_raster11[jj] = nodata_out
csmax_raster22[jj] = nodata_out
if (self.ssmname != ""):
ssm_raster[jj] = ssmLine[jj]
else:
ssm_raster[jj] = nodata_out
raster1a[jj] = normal[2]/(self.repeatTime/self.XSize/365.0/24.0/3600.0)*(normal[2]*yunit[1]-normal[1]*yunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2]));
raster1b[jj] = -normal[2]/(self.repeatTime/self.YSize/365.0/24.0/3600.0)*(normal[2]*xunit[1]-normal[1]*xunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2]));
raster2a[jj] = -normal[2]/(self.repeatTime/self.XSize/365.0/24.0/3600.0)*(normal[2]*yunit[0]-normal[0]*yunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2]));
raster2b[jj] = normal[2]/(self.repeatTime/self.YSize/365.0/24.0/3600.0)*(normal[2]*xunit[0]-normal[0]*xunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2]));
# pdb.set_trace()
poBand1.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster1.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32)
poBand2.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster2.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32)
if ((self.dhdxname != "")&(self.vxname != "")):
poBand1Off.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster11.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32)
poBand2Off.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster22.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32)
if ((self.dhdxname != "")&(self.srxname != "")):
poBand1Sch.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=sr_raster11.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32)
poBand2Sch.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=sr_raster22.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32)
if (self.csminxname != ""):
poBand1Min.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=csmin_raster11.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32)
poBand2Min.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=csmin_raster22.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32)
if (self.csmaxxname != ""):
poBand1Max.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=csmax_raster11.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32)
poBand2Max.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=csmax_raster22.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32)
if (self.ssmname != ""):
poBand1Msk.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=ssm_raster.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32)
if (self.dhdxname != ""):
poBand1RO2VX.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster1a.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64)
poBand2RO2VX.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster1b.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64)
poBand1RO2VY.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster2a.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64)
@ -654,17 +686,17 @@ class GeogridOptical():
poDstDS = None
if ((self.dhdxname != "")&(self.vxname != "")):
poDstDSOff = None
if ((self.dhdxname != "")&(self.srxname != "")):
poDstDSSch = None
if (self.csminxname != ""):
poDstDSMin = None
if (self.csmaxxname != ""):
poDstDSMax = None
if (self.ssmname != ""):
poDstDSMsk = None
if (self.dhdxname != ""):
poDstDSRO2VX = None
poDstDSRO2VY = None

1
contrib/geo_autoRIFT/geogrid/bindings/geogridmodule.cpp
View File

@ -403,6 +403,7 @@ PyObject* setRO2VYFilename(PyObject *self, PyObject *args)
return Py_BuildValue("i", 0);
}
PyObject* setLookSide(PyObject *self, PyObject *args)
{
uint64_t ptr;

307
contrib/geo_autoRIFT/geogrid/src/geogrid.cpp
View File

@ -60,24 +60,68 @@ void geoGrid::geogrid()
std::cout << "XLimits: " << xmin << " " << xmax << "\n";
std::cout << "YLimits: " << ymin << " " << ymax << "\n";
std::cout << "Extent in km: " << (xmax - xmin)/1000. << " " << (ymax - ymin)/1000. << "\n";
if (demname != "")
{
std::cout << "DEM: " << demname << "\n";
std::cout << "Velocities: " << vxname << " " << vyname << "\n";
std::cout << "Search Range: " << srxname << " " << sryname << "\n";
std::cout << "Chip Size Min: " << csminxname << " " << csminyname << "\n";
std::cout << "Chip Size Max: " << csmaxxname << " " << csmaxyname << "\n";
}
if (dhdxname != "")
{
std::cout << "Slopes: " << dhdxname << " " << dhdyname << "\n";
}
if (vxname != "")
{
std::cout << "Velocities: " << vxname << " " << vyname << "\n";
}
if (srxname != "")
{
std::cout << "Search Range: " << srxname << " " << sryname << "\n";
}
if (csminxname != "")
{
std::cout << "Chip Size Min: " << csminxname << " " << csminyname << "\n";
}
if (csmaxxname != "")
{
std::cout << "Chip Size Max: " << csmaxxname << " " << csmaxyname << "\n";
}
if (ssmname != "")
{
std::cout << "Stable Surface Mask: " << ssmname << "\n";
std::cout << "Output Nodata Value: " << nodata_out << "\n";
}
std::cout << "\nOutputs: \n";
std::cout << "Window locations: " << pixlinename << "\n";
if (dhdxname != "")
{
if (vxname != "")
{
std::cout << "Window offsets: " << offsetname << "\n";
std::cout << "Window search range: " << searchrangename << "\n";
std::cout << "Window chip size min: " << chipsizeminname << "\n";
std::cout << "Window chip size max: " << chipsizemaxname << "\n";
}
std::cout << "Window rdr_off2vel_x vector: " << ro2vx_name << "\n";
std::cout << "Window rdr_off2vel_y vector: " << ro2vy_name << "\n";
if (srxname != "")
{
std::cout << "Window search range: " << searchrangename << "\n";
}
}
if (csminxname != "")
{
std::cout << "Window chip size min: " << chipsizeminname << "\n";
}
if (csmaxxname != "")
{
std::cout << "Window chip size max: " << chipsizemaxname << "\n";
}
if (ssmname != "")
{
std::cout << "Window stable surface mask: " << stablesurfacemaskname << "\n";
}
std::cout << "Output Nodata Value: " << nodata_out << "\n";
std::cout << "\nStarting processing .... \n";
@ -223,6 +267,8 @@ void geoGrid::geogrid()
exit(101);
}
}
if (ssmname != "")
{
if (ssmDS == NULL)
{
std::cout << "Error opening stable surface mask file { " << ssmname << " }\n";
@ -230,6 +276,7 @@ void geoGrid::geogrid()
GDALDestroyDriverManager();
exit(101);
}
}
demDS->GetGeoTransform(geoTrans);
int demXSize = demDS->GetRasterXSize();
@ -251,7 +298,9 @@ void geoGrid::geogrid()
std::cout << "Ylimits : " << geoTrans[3] + (lOff + lCount) * geoTrans[5] << " "
<< geoTrans[3] + lOff * geoTrans[5] << "\n";
std::cout << "Dimensions of geogrid: " << pCount << " x " << lCount << "\n\n";
std::cout << "Origin index (in DEM) of geogrid: " << pOff << " " << lOff << "\n";
std::cout << "Dimensions of geogrid: " << pCount << " x " << lCount << "\n";
//Create GDAL Transformers
@ -330,6 +379,32 @@ void geoGrid::geogrid()
double raster2b[pCount];
// double raster2c[pCount];
GDALRasterBand *poBand1 = NULL;
GDALRasterBand *poBand2 = NULL;
GDALRasterBand *poBand1Off = NULL;
GDALRasterBand *poBand2Off = NULL;
GDALRasterBand *poBand1Sch = NULL;
GDALRasterBand *poBand2Sch = NULL;
GDALRasterBand *poBand1Min = NULL;
GDALRasterBand *poBand2Min = NULL;
GDALRasterBand *poBand1Max = NULL;
GDALRasterBand *poBand2Max = NULL;
GDALRasterBand *poBand1Msk = NULL;
GDALRasterBand *poBand1RO2VX = NULL;
GDALRasterBand *poBand1RO2VY = NULL;
GDALRasterBand *poBand2RO2VX = NULL;
GDALRasterBand *poBand2RO2VY = NULL;
GDALDataset *poDstDS = NULL;
GDALDataset *poDstDSOff = NULL;
GDALDataset *poDstDSSch = NULL;
GDALDataset *poDstDSMin = NULL;
GDALDataset *poDstDSMax = NULL;
GDALDataset *poDstDSMsk = NULL;
GDALDataset *poDstDSRO2VX = NULL;
GDALDataset *poDstDSRO2VY = NULL;
double nodata;
// double nodata_out;
if (vxname != "")
@ -342,7 +417,7 @@ void geoGrid::geogrid()
const char *pszFormat = "GTiff";
char **papszOptions = NULL;
std::string str = "";
double adfGeoTransform[6] = { xmin, (xmax-xmin)/pCount, 0, ymax, 0, (ymin-ymax)/lCount };
double adfGeoTransform[6] = { geoTrans[0] + pOff * geoTrans[1], geoTrans[1], 0, geoTrans[3] + lOff * geoTrans[5], 0, geoTrans[5]};
OGRSpatialReference oSRS;
char *pszSRS_WKT = NULL;
demSRS.exportToWkt( &pszSRS_WKT );
@ -353,7 +428,7 @@ void geoGrid::geogrid()
poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
if( poDriver == NULL )
exit(107);
GDALDataset *poDstDS;
// GDALDataset *poDstDS;
str = pixlinename;
const char * pszDstFilename = str.c_str();
@ -366,21 +441,22 @@ void geoGrid::geogrid()
// CPLFree( pszSRS_WKT );
GDALRasterBand *poBand1;
GDALRasterBand *poBand2;
// GDALRasterBand *poBand1;
// GDALRasterBand *poBand2;
poBand1 = poDstDS->GetRasterBand(1);
poBand2 = poDstDS->GetRasterBand(2);
poBand1->SetNoDataValue(nodata_out);
poBand2->SetNoDataValue(nodata_out);
if ((dhdxname != "")&(vxname != ""))
{
GDALDriver *poDriverOff;
poDriverOff = GetGDALDriverManager()->GetDriverByName(pszFormat);
if( poDriverOff == NULL )
exit(107);
GDALDataset *poDstDSOff;
// GDALDataset *poDstDSOff;
str = offsetname;
const char * pszDstFilenameOff = str.c_str();
@ -389,22 +465,25 @@ void geoGrid::geogrid()
poDstDSOff->SetGeoTransform( adfGeoTransform );
poDstDSOff->SetProjection( pszSRS_WKT );
// CPLFree( pszSRS_WKT );
// CPLFree( pszSRS_WKT );
GDALRasterBand *poBand1Off;
GDALRasterBand *poBand2Off;
// GDALRasterBand *poBand1Off;
// GDALRasterBand *poBand2Off;
poBand1Off = poDstDSOff->GetRasterBand(1);
poBand2Off = poDstDSOff->GetRasterBand(2);
poBand1Off->SetNoDataValue(nodata_out);
poBand2Off->SetNoDataValue(nodata_out);
}
if ((dhdxname != "")&(srxname != ""))
{
GDALDriver *poDriverSch;
poDriverSch = GetGDALDriverManager()->GetDriverByName(pszFormat);
if( poDriverSch == NULL )
exit(107);
GDALDataset *poDstDSSch;
// GDALDataset *poDstDSSch;
str = searchrangename;
const char * pszDstFilenameSch = str.c_str();
@ -415,20 +494,23 @@ void geoGrid::geogrid()
poDstDSSch->SetProjection( pszSRS_WKT );
// CPLFree( pszSRS_WKT );
GDALRasterBand *poBand1Sch;
GDALRasterBand *poBand2Sch;
// GDALRasterBand *poBand1Sch;
// GDALRasterBand *poBand2Sch;
poBand1Sch = poDstDSSch->GetRasterBand(1);
poBand2Sch = poDstDSSch->GetRasterBand(2);
poBand1Sch->SetNoDataValue(nodata_out);
poBand2Sch->SetNoDataValue(nodata_out);
}
if (csminxname != "")
{
GDALDriver *poDriverMin;
poDriverMin = GetGDALDriverManager()->GetDriverByName(pszFormat);
if( poDriverMin == NULL )
exit(107);
GDALDataset *poDstDSMin;
// GDALDataset *poDstDSMin;
str = chipsizeminname;
const char * pszDstFilenameMin = str.c_str();
@ -439,21 +521,24 @@ void geoGrid::geogrid()
poDstDSMin->SetProjection( pszSRS_WKT );
// CPLFree( pszSRS_WKT );
GDALRasterBand *poBand1Min;
GDALRasterBand *poBand2Min;
// GDALRasterBand *poBand1Min;
// GDALRasterBand *poBand2Min;
poBand1Min = poDstDSMin->GetRasterBand(1);
poBand2Min = poDstDSMin->GetRasterBand(2);
poBand1Min->SetNoDataValue(nodata_out);
poBand2Min->SetNoDataValue(nodata_out);
}
if (csmaxxname != "")
{
GDALDriver *poDriverMax;
poDriverMax = GetGDALDriverManager()->GetDriverByName(pszFormat);
if( poDriverMax == NULL )
exit(107);
GDALDataset *poDstDSMax;
// GDALDataset *poDstDSMax;
str = chipsizemaxname;
const char * pszDstFilenameMax = str.c_str();
@ -464,22 +549,25 @@ void geoGrid::geogrid()
poDstDSMax->SetProjection( pszSRS_WKT );
// CPLFree( pszSRS_WKT );
GDALRasterBand *poBand1Max;
GDALRasterBand *poBand2Max;
// GDALRasterBand *poBand1Max;
// GDALRasterBand *poBand2Max;
poBand1Max = poDstDSMax->GetRasterBand(1);
poBand2Max = poDstDSMax->GetRasterBand(2);
poBand1Max->SetNoDataValue(nodata_out);
poBand2Max->SetNoDataValue(nodata_out);
}
if (ssmname != "")
{
GDALDriver *poDriverMsk;
poDriverMsk = GetGDALDriverManager()->GetDriverByName(pszFormat);
if( poDriverMsk == NULL )
exit(107);
GDALDataset *poDstDSMsk;
// GDALDataset *poDstDSMsk;
str = stablesurfacemaskname;
const char * pszDstFilenameMsk = str.c_str();
@ -490,18 +578,21 @@ void geoGrid::geogrid()
poDstDSMsk->SetProjection( pszSRS_WKT );
// CPLFree( pszSRS_WKT );
GDALRasterBand *poBand1Msk;
// GDALRasterBand *poBand1Msk;
poBand1Msk = poDstDSMsk->GetRasterBand(1);
poBand1Msk->SetNoDataValue(nodata_out);
}
if (dhdxname != "")
{
GDALDriver *poDriverRO2VX;
poDriverRO2VX = GetGDALDriverManager()->GetDriverByName(pszFormat);
if( poDriverRO2VX == NULL )
exit(107);
GDALDataset *poDstDSRO2VX;
// GDALDataset *poDstDSRO2VX;
str = ro2vx_name;
const char * pszDstFilenameRO2VX = str.c_str();
@ -510,25 +601,24 @@ void geoGrid::geogrid()
poDstDSRO2VX->SetGeoTransform( adfGeoTransform );
poDstDSRO2VX->SetProjection( pszSRS_WKT );
// CPLFree( pszSRS_WKT );
// CPLFree( pszSRS_WKT );
GDALRasterBand *poBand1RO2VX;
GDALRasterBand *poBand2RO2VX;
// GDALRasterBand *poBand3Los;
// GDALRasterBand *poBand1RO2VX;
// GDALRasterBand *poBand2RO2VX;
// GDALRasterBand *poBand3Los;
poBand1RO2VX = poDstDSRO2VX->GetRasterBand(1);
poBand2RO2VX = poDstDSRO2VX->GetRasterBand(2);
// poBand3Los = poDstDSLos->GetRasterBand(3);
// poBand3Los = poDstDSLos->GetRasterBand(3);
poBand1RO2VX->SetNoDataValue(nodata_out);
poBand2RO2VX->SetNoDataValue(nodata_out);
// poBand3Los->SetNoDataValue(nodata_out);
// poBand3Los->SetNoDataValue(nodata_out);
GDALDriver *poDriverRO2VY;
poDriverRO2VY = GetGDALDriverManager()->GetDriverByName(pszFormat);
if( poDriverRO2VY == NULL )
exit(107);
GDALDataset *poDstDSRO2VY;
// GDALDataset *poDstDSRO2VY;
str = ro2vy_name;
const char * pszDstFilenameRO2VY = str.c_str();
@ -537,17 +627,24 @@ void geoGrid::geogrid()
poDstDSRO2VY->SetGeoTransform( adfGeoTransform );
poDstDSRO2VY->SetProjection( pszSRS_WKT );
CPLFree( pszSRS_WKT );
// CPLFree( pszSRS_WKT );
GDALRasterBand *poBand1RO2VY;
GDALRasterBand *poBand2RO2VY;
// GDALRasterBand *poBand3Alt;
// GDALRasterBand *poBand1RO2VY;
// GDALRasterBand *poBand2RO2VY;
// GDALRasterBand *poBand3Alt;
poBand1RO2VY = poDstDSRO2VY->GetRasterBand(1);
poBand2RO2VY = poDstDSRO2VY->GetRasterBand(2);
// poBand3Alt = poDstDSAlt->GetRasterBand(3);
// poBand3Alt = poDstDSAlt->GetRasterBand(3);
poBand1RO2VY->SetNoDataValue(nodata_out);
poBand2RO2VY->SetNoDataValue(nodata_out);
// poBand3Alt->SetNoDataValue(nodata_out);
// poBand3Alt->SetNoDataValue(nodata_out);
}
CPLFree( pszSRS_WKT );
// ground range and azimuth pixel size
@ -865,6 +962,8 @@ void geoGrid::geogrid()
double los[3];
double alt[3];
double normal[3];
double cross[3];
double cross_check;
double dopfact;
double height;
@ -1110,18 +1209,19 @@ void geoGrid::geogrid()
raster2a[jj] = nodata_out;
raster2b[jj] = nodata_out;
// raster2c[jj] = nodata_out;
}
else
{
raster1[jj] = rgind;
raster2[jj] = azind;
if (vxname == "")
if (dhdxname != "")
{
raster11[jj] = 0.;
raster22[jj] = 0.;
}
else if (vel[0] == nodata)
if (vxname != "")
{
if (vel[0] == nodata)
{
raster11[jj] = 0.;
raster22[jj] = 0.;
@ -1132,12 +1232,30 @@ void geoGrid::geogrid()
raster22[jj] = std::round(dot_C(vel,temp)*dt/norm_C(alt)/365.0/24.0/3600.0*1);
}
if (srxname == "")
{
sr_raster11[jj] = nodata_out;
sr_raster22[jj] = nodata_out;
}
else if (vel[0] == nodata)
cross_C(los,temp,cross);
unitvec_C(cross, cross);
cross_check = std::abs(std::acos(dot_C(normal,cross))/deg2rad-90.0);
if (cross_check > 1.0)
{
raster1a[jj] = normal[2]/(dt/dr/365.0/24.0/3600.0)*(normal[2]*temp[1]-normal[1]*temp[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2]));
raster1b[jj] = -normal[2]/(dt/norm_C(alt)/365.0/24.0/3600.0)*(normal[2]*los[1]-normal[1]*los[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2]));
raster2a[jj] = -normal[2]/(dt/dr/365.0/24.0/3600.0)*(normal[2]*temp[0]-normal[0]*temp[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2]));
raster2b[jj] = normal[2]/(dt/norm_C(alt)/365.0/24.0/3600.0)*(normal[2]*los[0]-normal[0]*los[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2]));
}
else
{
raster1a[jj] = nodata_out;
raster1b[jj] = nodata_out;
raster2a[jj] = nodata_out;
raster2b[jj] = nodata_out;
}
if (srxname != "")
{
if ((vxname != "")&(vel[0] == nodata))
{
sr_raster11[jj] = 0;
sr_raster22[jj] = 0;
@ -1163,39 +1281,32 @@ void geoGrid::geogrid()
sr_raster22[jj] = 1;
}
}
}
}
if (csminxname != "")
{
csmin_raster11[jj] = csminxLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_grd;
csmin_raster22[jj] = csminyLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_azm;
}
else
{
csmin_raster11[jj] = nodata_out;
csmin_raster22[jj] = nodata_out;
}
if (csmaxxname != "")
{
csmax_raster11[jj] = csmaxxLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_grd;
csmax_raster22[jj] = csmaxyLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_azm;
}
else
{
csmax_raster11[jj] = nodata_out;
csmax_raster22[jj] = nodata_out;
}
if (ssmname != "")
{
ssm_raster[jj] = ssmLine[jj];
}
else
{
ssm_raster[jj] = nodata_out;
}
raster1a[jj] = normal[2]/(dt/dr/365.0/24.0/3600.0)*(normal[2]*temp[1]-normal[1]*temp[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2]));
raster1b[jj] = -normal[2]/(dt/norm_C(alt)/365.0/24.0/3600.0)*(normal[2]*los[1]-normal[1]*los[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2]));
raster2a[jj] = -normal[2]/(dt/dr/365.0/24.0/3600.0)*(normal[2]*temp[0]-normal[0]*temp[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2]));
raster2b[jj] = normal[2]/(dt/norm_C(alt)/365.0/24.0/3600.0)*(normal[2]*los[0]-normal[0]*los[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2]));
// raster1a[jj] = los[0]*dt/dr/365.0/24.0/3600.0;
// raster1b[jj] = los[1]*dt/dr/365.0/24.0/3600.0;
@ -1212,68 +1323,112 @@ void geoGrid::geogrid()
// std::cout << raster1[ii][jj] << "\n";
}
poBand1->RasterIO( GF_Write, 0, ii, pCount, 1,
raster1, pCount, 1, GDT_Int32, 0, 0 );
poBand2->RasterIO( GF_Write, 0, ii, pCount, 1,
raster2, pCount, 1, GDT_Int32, 0, 0 );
if ((dhdxname != "")&(vxname != ""))
{
poBand1Off->RasterIO( GF_Write, 0, ii, pCount, 1,
raster11, pCount, 1, GDT_Int32, 0, 0 );
poBand2Off->RasterIO( GF_Write, 0, ii, pCount, 1,
raster22, pCount, 1, GDT_Int32, 0, 0 );
}
if ((dhdxname != "")&(srxname != ""))
{
poBand1Sch->RasterIO( GF_Write, 0, ii, pCount, 1,
sr_raster11, pCount, 1, GDT_Int32, 0, 0 );
poBand2Sch->RasterIO( GF_Write, 0, ii, pCount, 1,
sr_raster22, pCount, 1, GDT_Int32, 0, 0 );
}
if (csminxname != "")
{
poBand1Min->RasterIO( GF_Write, 0, ii, pCount, 1,
csmin_raster11, pCount, 1, GDT_Int32, 0, 0 );
poBand2Min->RasterIO( GF_Write, 0, ii, pCount, 1,
csmin_raster22, pCount, 1, GDT_Int32, 0, 0 );
}
if (csmaxxname != "")
{
poBand1Max->RasterIO( GF_Write, 0, ii, pCount, 1,
csmax_raster11, pCount, 1, GDT_Int32, 0, 0 );
poBand2Max->RasterIO( GF_Write, 0, ii, pCount, 1,
csmax_raster22, pCount, 1, GDT_Int32, 0, 0 );
}
if (ssmname != "")
{
poBand1Msk->RasterIO( GF_Write, 0, ii, pCount, 1,
ssm_raster, pCount, 1, GDT_Int32, 0, 0 );
}
if (dhdxname != "")
{
poBand1RO2VX->RasterIO( GF_Write, 0, ii, pCount, 1,
raster1a, pCount, 1, GDT_Float64, 0, 0 );
poBand2RO2VX->RasterIO( GF_Write, 0, ii, pCount, 1,
raster1b, pCount, 1, GDT_Float64, 0, 0 );
// poBand3Los->RasterIO( GF_Write, 0, ii, pCount, 1,
// raster1c, pCount, 1, GDT_Float64, 0, 0 );
// poBand3Los->RasterIO( GF_Write, 0, ii, pCount, 1,
// raster1c, pCount, 1, GDT_Float64, 0, 0 );
poBand1RO2VY->RasterIO( GF_Write, 0, ii, pCount, 1,
raster2a, pCount, 1, GDT_Float64, 0, 0 );
poBand2RO2VY->RasterIO( GF_Write, 0, ii, pCount, 1,
raster2b, pCount, 1, GDT_Float64, 0, 0 );
// poBand3Alt->RasterIO( GF_Write, 0, ii, pCount, 1,
// raster2c, pCount, 1, GDT_Float64, 0, 0 );
// poBand3Alt->RasterIO( GF_Write, 0, ii, pCount, 1,
// raster2c, pCount, 1, GDT_Float64, 0, 0 );
}
}
/* Once we're done, close properly the dataset */
GDALClose( (GDALDatasetH) poDstDS );
if ((dhdxname != "")&(vxname != ""))
{
/* Once we're done, close properly the dataset */
GDALClose( (GDALDatasetH) poDstDSOff );
}
if ((dhdxname != "")&(srxname != ""))
{
/* Once we're done, close properly the dataset */
GDALClose( (GDALDatasetH) poDstDSSch );
}
if (csminxname != "")
{
/* Once we're done, close properly the dataset */
GDALClose( (GDALDatasetH) poDstDSMin );
}
if (csmaxxname != "")
{
/* Once we're done, close properly the dataset */
GDALClose( (GDALDatasetH) poDstDSMax );
}
if (ssmname != "")
{
/* Once we're done, close properly the dataset */
GDALClose( (GDALDatasetH) poDstDSMsk );
}
if (dhdxname != "")
{
/* Once we're done, close properly the dataset */
GDALClose( (GDALDatasetH) poDstDSRO2VX );
/* Once we're done, close properly the dataset */
GDALClose( (GDALDatasetH) poDstDSRO2VY );
}
GDALClose(demDS);