注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

阿弥陀佛

街树飘影未见尘 潭月潜水了无声 般若观照心空静...

 
 
 

日志

 
 
关于我

一直从事气象预报、服务建模实践应用。 注重气象物理场、实况场、地理信息、本体知识库、分布式气象内容管理系统建立。 对Barnes客观分析, 小波,计算神经网络、信任传播、贝叶斯推理、专家系统、网络本体语言有一定体会。 一直使用Java、Delphi、Prolog、SQL编程。

网易考拉推荐

使用雙線性插值完成Lambert到LatLon投影的轉換  

2014-12-02 19:58:30|  分类: Scala |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
package interpolation
import java.io.{ File, BufferedReader, ByteArrayInputStream, InputStreamReader }
import ucar.nc2.ncml.{ NcMLReader }
import ucar.nc2.{ NetcdfFile, NCdumpW }
import scala.xml._
import java.nio.charset.{ Charset, CodingErrorAction }
object LapsNcml {
  def main(args: Array[String]) {
    makeNcFileFromNcml
  }
  def makeNcFileFromNcml {
    val ncFnm = "/wk/laps/lapstest1.nc"
    val ncml = ncmlTemp
    val inputStream = new ByteArrayInputStream(ncml.getBytes("UTF-8"))
    val decoder = Charset.forName("UTF-8").newDecoder();
    decoder.onMalformedInput(CodingErrorAction.IGNORE);
    val char_input = new InputStreamReader(inputStream, decoder)
    val ncd = NcMLReader.readNcML(char_input, null)
    val ncdnew = ucar.nc2.FileWriter.writeToFile(ncd, ncFnm, true)
    ncd.close
    ncdnew.close
  }
  def ncmlTemp = {
    "<?xml version=\"1.0\" encoding=\"UTF-8\"?>" +
      <netcdf xmlns="http://www.unidata.ucar.edu/namespaces/netcdf/ncml-2.2">
        <dimension name="time" length="1" isUnlimited="true"/>
        <dimension name="lat" length="517"/>
        <dimension name="lon" length="947"/>
        <attribute name="title" type="String" value={ "LAPS 地面气象要素格点" }/>
        <attribute name="missing_value" value="NaN" type="float"/>
        <variable name="lat" shape="lat" type="float">
          <attribute name="units" type="String" value="degrees_north"/>
          <values start="10.160842895507812" increment="0.1"/>
        </variable>
        <variable name="lon" shape="lon" type="float">
          <attribute name="units" type="String" value="degrees_east"/>
          <values start="54.624755859375" increment="0.1"/>
        </variable>
        <variable name="time" shape="time" type="double">
          <attribute name="long_name" type="String"/>
          <attribute name="units" type="String" value="seconds since 0001-01-01 00:00:00"/>
        </variable>
        <variable name="t" shape="time lat lon" type="float">
          <attribute name="navigation_dim" value="nav"/>
          <attribute name="record" value="valtime, reftime"/>
          <attribute name="_FillValue" type="float" value="1.0E37"/>
          <attribute name="long_name" value="氣溫"/>
          <attribute name="units" value="kelvin"/>
          <attribute name="valid_range" type="float" value="0.0 350.0"/>
          <attribute name="LAPS_var" value="T"/>
          <attribute name="lvl_coord" value="AGL "/>
          <attribute name="LAPS_units" value="K"/>
        </variable>
        <variable name="rh" shape="time lat lon" type="float">
          <attribute name="navigation_dim" value="nav"/>
          <attribute name="record" value="valtime, reftime"/>
          <attribute name="_FillValue" type="float" value="1.0E37"/>
          <attribute name="long_name" value="相對溼度"/>
          <attribute name="units" value="meter"/>
          <attribute name="valid_range" type="float" value="-20000.0 20000.0"/>
          <attribute name="LAPS_var" value="RH"/>
          <attribute name="lvl_coord" value="AGL"/>
          <attribute name="LAPS_units" value="M"/>
        </variable>
        <variable name="ps" shape="time lat lon" type="float">
          <attribute name="navigation_dim" value="nav"/>
          <attribute name="record" value="valtime, reftime"/>
          <attribute name="_FillValue" type="float" value="1.0E37"/>
          <attribute name="long_name" value="本站氣壓"/>
          <attribute name="units" value="pascal"/>
          <attribute name="valid_range" type="float" value="0.0 120000.0"/>
          <attribute name="LAPS_var" value="PS"/>
          <attribute name="lvl_coord" value="AGL"/>
          <attribute name="LAPS_units" value="PA"/>
        </variable>
      </netcdf>.toString
  }
}
===========================================================================
package interpolation
import com.github.nscala_time.time.Imports._
import java.io.{ File, BufferedReader, ByteArrayInputStream, InputStreamReader }
import ucar.nc2.ncml.{ NcMLReader }
import ucar.nc2.{ NetcdfFile, NCdumpW, Variable, NetcdfFileWriter }
import java.nio.charset.{ Charset, CodingErrorAction }
import ucar.unidata.geoloc.projection.LambertConformal
import ucar.unidata.geoloc.ProjectionPointImpl
object Lambert2LatLon {
  val nf = NetcdfFile.open("/home/hxf/static.nest7grid")
  val Ny = nf.readSection("Ny").toString.trim.toInt
  val Nx = nf.readSection("Nx").toString.trim.toInt
  val latmat = getMat(nf, "lat", Ny, Nx)
  val lonmat = getMat(nf, "lon", Ny, Nx)
  nf.close

  val nRow = 517
  val nCol = 947
  val step = 0.1
  //LatLon起始經緯度
  val lat0 = 10.160842895507812
  val lon0 = 54.624755859375
  def main(args: Array[String]) {
    val fnmSrc = "/wk/laps/143292300.lsx"
    val fnmDst = "/wk/laps/lapstest1.nc"
    val time = "2014-03-29"
    interpolation(fnmSrc, fnmDst, time)
  }

  def interpolation(src: String, dst: String, time: String) {
    val nfsrc = NetcdfFile.open(src)
    val nfdst = NetcdfFileWriter.openExisting(dst)
    val dt = new DateTime(time)
    appendTime(dt, nfdst)
    interEle(nfsrc, nfdst, "t")
    interEle(nfsrc, nfdst, "ps")
    interEle(nfsrc, nfdst, "rh")
    nfdst.flush()
    nfdst.close()
  }
  type matType = Array[Array[Float]]

  def getMat(nf: NetcdfFile, nm: String, ny: Int, nx: Int): matType = {
    val a = nf.readSection(nm).copyTo1DJavaArray().asInstanceOf[Array[Float]]
    val mat = Array.fill(ny, nx)(0.0f)
    var k = 0
    for (i <- 0 to ny - 1) {
      Array.copy(a, k, mat(i), 0, nx)
      k += nx
    }
    mat
  }
  def appendTime(time0: DateTime, nc: NetcdfFileWriter) {
    val tvar = time0.getMillis
    val t0 = new DateTime(1, 1, 1, 0, 0, 0, 0).getMillis
    val dt = (tvar - t0 + 1000.0) / 1000.0
    val variable = nc.findVariable("time")
    val time_val = Array(dt)
    val a = ucar.ma2.Array.factory(time_val)
    nc.write(variable, a)
  }

  def interEle(nfsrc: NetcdfFile, nfdst: NetcdfFileWriter, eleNm: String) {
    val elesrc = getMat(nfsrc, eleNm, Ny, Nx)
    val eledst = getMat(nfdst.getNetcdfFile, eleNm, nRow, nCol)
    for (
      y <- 0 to nRow - 1;
      x <- 0 to nCol - 1;
      latDst = y * step + lat0;
      lonDst = x * step + lon0;
      ij = Lapslatlon2ij(latDst, lonDst);
      i = ij._1;
      j = ij._2;
      i1 = i + 1;
      j1 = j + 1;
      if (i >= 0 && j >= 0 && i1 < Ny && j1 < Nx);
      x1 = lonmat(i)(j);
      y1 = latmat(i)(j);
      x2 = lonmat(i1)(j1);
      y2 = latmat(i1)(j1);
      q11 = elesrc(i)(j);
      q21 = elesrc(i)(j1);
      q12 = elesrc(i1)(j);
      q22 = elesrc(i1)(j1);
      bi = bilinear(latDst, lonDst, x1, x2, y1, y2, q11, q12, q21, q22)
    ) {
      eledst(y)(x) = bi
    }
    val a = ucar.ma2.Array.factory(eledst)
    val v = nfdst.findVariable(eleNm)
    val b = a.reshape(v.getShape)
    nfdst.write(v, b)
  }

  def Lapslatlon2ij(lat: Double, lon: Double) = {
    val Dxy = 3
    //中國Lambert投影格點場起始經緯度, 注意起始經度 > 其中的lon矩陣中的最小值54.624755859375
    val La1 = 10.160843
    val Lo1 = 78.986084
    val lambert = new LambertConformal(37.0, 102.0, 30.0, 60.0) //中國Lambert投影參數
    val xy0 = lambert.latLonToProj(La1, Lo1)
    val xy = lambert.latLonToProj(lat, lon)
    val x = -xy0.getX + xy.getX
    val y = -xy0.getY + xy.getY
    val i = (y / Dxy).toInt
    val j = (x / Dxy).toInt
    (i, j)
  }

  /**
   * 雙線性插值. 見http://en.wikipedia.org/wiki/Bilinear_interpolation
   */
  def bilinear(lat: Double, lon: Double,
    x1: Double, x2: Double, y1: Double, y2: Double,
    q11: Double, q12: Double, q21: Double, q22: Double) = {
    val d = 1 / ((x2 - x1) * (y2 - y1))
    val dq11 = q11 * (x2 - lon) * (y2 - lat)
    val dq21 = q21 * (lon - x1) * (y2 - lat)
    val dq12 = q12 * (x2 - lon) * (lat - y1)
    val dq22 = q22 * (lon - x1) * (lat - y1)
    val v = (d * (dq11 + dq21 + dq12 + dq22)).toFloat
    v
  }
}
使用雙線性插值完成Lambert到LatLon投影的轉換 - 险峰 - 阿弥陀佛
  评论这张
 
阅读(448)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017