A bug has been fixed: there was a confusion between the mean adiabatic lapse...

A bug has been fixed: there was a confusion between the mean adiabatic lapse rate and the dry one. Small code cleanup in DEMObject (with a better toString() method), ARPSIO and ProcShade.

meteoio/dataClasses/DEMObject.cc

@@ -679,7 +679,7 @@ double DEMObject::getHorizon(const Coords& point, const double& bearing) const
 		if (new_altitude==mio::IOUtils::nodata) break; //we stop at nodata cells
 
 		const double DeltaH = new_altitude - cell_alt;
-		const double distance = sqrt( (double)( (ix2-ix1)*(ix2-ix1) + (iy2-iy1)*(iy2-iy1)) ) * cellsize;
+		const double distance = sqrt( (double)( Optim::pow2(ix2-ix1) + Optim::pow2(iy2-iy1)) ) * cellsize;
 		const double tan_angle = DeltaH/distance;
 		if (tan_angle>horizon_tan_angle) horizon_tan_angle = tan_angle;
 
@@ -1059,6 +1059,20 @@ double DEMObject::safeGet(const int i, const int j)
 	return grid2D((unsigned)i, (unsigned)j);
 }
 
+const std::string DEMObject::toString() const {
+	std::ostringstream os;
+	os << "<DEMObject>\n";
+	os << llcorner.toString();
+	os << grid2D.getNx() << " x " << grid2D.getNy() << " @ " << cellsize << "m\t[ " << min_altitude << " - " << max_altitude << " ]\n";
+	//os << grid2D.toString();
+	os << "Slope: " << slope.getNx() << " x " << slope.getNy() << "\t[ " << min_slope << " - " << max_slope << " ]\n";
+	os << "Azi: " << azi.getNx() << " x " << azi.getNy() << "\n";
+	os << "Curvature: " << curvature.getNx() << " x " << curvature.getNy()<< "\t[ " << min_curvature << " - " << max_curvature << " ]\n";
+	os << "Nx: " << Nx.getNx() << " x " << Nx.getNy() << " , Ny: " << Ny.getNx() << " x " << Ny.getNy() << " , Nz: " << Nz.getNx() << " x " << Nz.getNy() << "\n";
+	os << "</DEMObject>\n";
+	return os.str();
+}
+
 std::iostream& operator<<(std::iostream& os, const DEMObject& dem) {
 	operator<<(os, *((Grid2DObject*)&dem));
 

meteoio/dataClasses/DEMObject.h

@@ -129,6 +129,7 @@ class DEMObject : public Grid2DObject {
 		bool operator==(const DEMObject& in) const; ///<Operator that tests for equality
 		bool operator!=(const DEMObject& in) const; ///<Operator that tests for inequality
 
+		const std::string toString() const;
 		friend std::iostream& operator<<(std::iostream& os, const DEMObject& dem);
 		friend std::iostream& operator>>(std::iostream& is, DEMObject& dem);
 

meteoio/meteoFilters/ProcShade.cc

@@ -52,7 +52,7 @@ void ProcShade::process(const unsigned int& param, const std::vector<MeteoData>&
                         std::vector<MeteoData>& ovec)
 {
 	ovec = ivec;
-	if (ovec.size()==0) return;
+	if (ovec.empty()) return;
 	
 	const string stationHash = ovec[0].meta.getHash();
 	
@@ -71,13 +71,8 @@ void ProcShade::process(const unsigned int& param, const std::vector<MeteoData>&
 		//now look for our specific station hash
 		mask = masks.find( stationHash );
 		if (mask==masks.end()) {
-			Coords position( ovec[0].meta.position );
-			if (!dem.gridify(position)) {
-				const string msg = "In filter '"+block_name+"', station '"+ovec[0].meta.stationID+"' "+position.toString(Coords::LATLON)+" is not included in the DEM "+dem.llcorner.toString(Coords::LATLON);
-				throw NoDataException(msg, AT);
-			}
 			std::vector< std::pair<double,double> > tmp_mask;
-			dem.getHorizon(position, 10., tmp_mask);
+			computeMask(dem, ovec[0].meta, tmp_mask);
 			masks[ stationHash ] = tmp_mask;
 			mask = masks.find( stationHash);
 		}
@@ -207,11 +202,30 @@ void ProcShade::readMask(const std::string& filter, const std::string& filename,
 	std::sort(o_mask.begin(), o_mask.end(), sort_pred());
 }
 
+void ProcShade::computeMask(const DEMObject& i_dem, const StationData& sd, std::vector< std::pair<double,double> > &o_mask, const bool& dump_mask)
+{
+	o_mask.clear();
+	Coords position( sd.position );
+	if (!i_dem.gridify(position)) {
+		const string msg = "In filter 'SHADE', station '"+sd.stationID+"' "+position.toString(Coords::LATLON)+" is not included in the DEM "+i_dem.llcorner.toString(Coords::LATLON);
+		throw NoDataException(msg, AT);
+	}
+	
+	i_dem.getHorizon(position, 10., o_mask); //by 10deg increments
+	if (dump_mask) {
+		std::cout << "Horizon mask for station '" << sd.stationID << "'\n";
+		for (size_t ii=0; ii<o_mask.size(); ii++)
+			std::cout << o_mask[ii].first << " " << o_mask[ii].second << "\n";
+		std::cout << "\n";
+	}
+}
+
 void ProcShade::parse_args(const std::vector<std::string>& vec_args)
 {
 	const size_t nrArgs = vec_args.size();
 	if (nrArgs==0) { //compute from DEM
 		IOHandler io(cfg);
+		dem.setUpdatePpt( DEMObject::NO_UPDATE ); //we only need the elevations
 		io.readDEM(dem);
 	} else if (nrArgs==1) {
 		const std::string root_path( cfg.getConfigRootDir() );

meteoio/meteoFilters/ProcShade.h

@@ -62,10 +62,11 @@ class ProcShade : public ProcessingBlock {
 
 		virtual void process(const unsigned int& param, const std::vector<MeteoData>& ivec,
 		                     std::vector<MeteoData>& ovec);
+		
+		static void computeMask(const DEMObject& i_dem, const StationData& sd, std::vector< std::pair<double,double> > &o_mask, const bool& dump_mask=false);
 
 	private:
 		static void readMask(const std::string& filter, const std::string& filename, std::vector< std::pair<double,double> > &o_mask);
-		
 		void parse_args(const std::vector<std::string>& vec_args);
 		double getMaskElevation(const std::vector< std::pair<double,double> > &mask, const double& azimuth) const;
 

meteoio/meteoLaws/Atmosphere.cc

@@ -63,8 +63,8 @@ double Atmosphere::blkBody_Radiation(const double& ea, const double& T) {
 * @return standard pressure (Pa)
 */
 double Atmosphere::stdAirPressure(const double& altitude) {
-	const double expo = Cst::gravity / (Cst::dry_adiabatique_lapse_rate * Cst::gaz_constant_dry_air);
-	const double p = Cst::std_press * pow( 1. - ( (Cst::dry_adiabatique_lapse_rate * Cst::earth_R0 * altitude) / (Cst::std_temp * (Cst::earth_R0 + altitude)) ), expo );
+	const double expo = Cst::gravity / (Cst::mean_adiabatique_lapse_rate * Cst::gaz_constant_dry_air);
+	const double p = Cst::std_press * pow( 1. - ( (Cst::mean_adiabatique_lapse_rate * Cst::earth_R0 * altitude) / (Cst::std_temp * (Cst::earth_R0 + altitude)) ), expo );
 	return p;
 }
 

meteoio/meteoLaws/Meteoconst.h

@@ -38,7 +38,7 @@ namespace Cst {
 	const double gravity = 9.80665; // (m s-2)
 	const double std_press = 101325.; // (Pa) at sea level
 	const double std_temp = 288.15; // (K) at sea level
-	const double dry_adiabatique_lapse_rate = 0.0065; // (K/m)
+	const double mean_adiabatique_lapse_rate = 0.0065; // (K/m)
 
 	const double gaz_constant_dry_air = 287.058; // (J kg-1 K-1)
 	const double gaz_constant_water_vapor = 461.9; // (J kg-1 K-1)

meteoio/plugins/ARPSIO.cc

@@ -542,9 +542,12 @@ void ARPSIO::readGridLayer(FILE* &fin, const std::string& filename, const std::s
 				grid(ix,iy) = tmp;
 			} else {
 				char dummy[ARPS_MAX_LINE_LENGTH];
-				(void)fscanf(fin,"%s",dummy);
+				const int status = fscanf(fin,"%s",dummy);
 				fclose(fin);
-				throw InvalidFormatException("Fail to read data layer for parameter '"+parameter+"' in file '"+filename+"', instead read: '"+string(dummy)+"'", AT);
+				string msg = "Fail to read data layer for parameter '"+parameter+"' in file '"+filename+"'";
+				if (status>0) msg += ", instead read: '"+string(dummy)+"'";
+				
+				throw InvalidFormatException(msg, AT);
 			}
 		}
 	}
@@ -553,24 +556,19 @@ void ARPSIO::readGridLayer(FILE* &fin, const std::string& filename, const std::s
 void ARPSIO::moveToMarker(FILE* &fin, const std::string& filename, const std::string& marker)
 {
 	char dummy[ARPS_MAX_LINE_LENGTH];
-	int nb_elems = 0;
 	do {
-		nb_elems = fscanf(fin," %[^\t\n] ",dummy); //HACK: possible buffer overflow
+		const int nb_elems = fscanf(fin," %[^\t\n] ",dummy); //HACK: possible buffer overflow
+		if (nb_elems==0) {
+			fclose(fin);
+			throw InvalidFormatException("Matching failure in file "+filename+" when looking for "+marker, AT);
+		}
 		if (dummy[0]>='A' && dummy[0]<='z') { //this is a "marker" line
 			if (strncmp(dummy,marker.c_str(), marker.size()) == 0) return;
 		}
-	} while (!feof(fin) && nb_elems!=0);
+	} while (!feof(fin));
 	
-	if (feof(fin)) {
-		fclose(fin);
-		const std::string message = "End of file "+filename+" should NOT have been reached when looking for "+marker;
-		throw InvalidFormatException(message, AT);
-	}
-	if (nb_elems==0) { //there was a line that did not match the format
-		fclose(fin);
-		const std::string message = "Matching failure in file "+filename+" when looking for "+marker;
-		throw InvalidFormatException(message, AT);
-	}
+	fclose(fin);
+	throw InvalidFormatException("End of file "+filename+" should NOT have been reached when looking for "+marker, AT);
 }
 
 } //namespace

meteoio/plugins/NetCDFIO.cc

@@ -439,11 +439,11 @@ void NetCDFIO::readDEM(DEMObject& dem_out)
 			read2DGrid_internal(p, filename, MeteoGrids::P) && 
 			read2DGrid_internal(ta, filename, MeteoGrids::TA)) {
 				dem_out.set(p, IOUtils::nodata);
-				const double k = Cst::gravity / (Cst::dry_adiabatique_lapse_rate * Cst::gaz_constant_dry_air);
+				const double k = Cst::gravity / (Cst::mean_adiabatique_lapse_rate * Cst::gaz_constant_dry_air);
 				const double k_inv = 1./k;
 				for (size_t ii=0; ii<(dem_out.getNx()*dem_out.getNy()); ii++) {
 					const double K = pow(p(ii)/p_sea(ii), k_inv);
-					dem_out(ii) = ta(ii)*Cst::earth_R0*(1.-K) / (Cst::dry_adiabatique_lapse_rate * Cst::earth_R0 - ta(ii)*(1.-K));
+					dem_out(ii) = ta(ii)*Cst::earth_R0*(1.-K) / (Cst::mean_adiabatique_lapse_rate * Cst::earth_R0 - ta(ii)*(1.-K));
 				}
 				return;
 			}