ラベル GIS の投稿を表示しています。 すべての投稿を表示
ラベル GIS の投稿を表示しています。 すべての投稿を表示

2026年3月30日月曜日

boost::geometry::extensions が boost_1_90_0 で工事中

boost::geometry の 拡張アルゴリズムを利用しておりまして、boost がリリースされる毎に boost::geometry のリポジトリから対応するものを持ってきて使っておりました。
no_rescale_policy というのが廃止されるとの事で、dissolve のルーチンも graph ベースのルーチンに絶賛改良の工事中でコンパイルが通りませんでした。
ちょちょっと、前のバージョンのものを流用すれば良いと高を括っていたら、結構な変更が入っていて、全然コンパイルが通らない…

という事で、一時的にでも手直しをして利用する事にしました。
一応 dissolve のテストは全部通したいと思って GPT-5.3-Codex にお手伝いいただいて、修正をかけようとしたら、
GPTが以下のプロセスに入って終了しない。

  •  修正アルゴリズムを考える
  •  修正を試みる 
  • テストをする 
  • テストの成績が平行または悪化
  • 修正を巻き戻す
GPTはダメだと思って Opus-4.6 にバトンタッチして修正をさせる事に…
通らないテストについて聞いてみると、自己交差をした図形が自動修復されるかテストしているとの事
正解をsvgで図形化して見せてもらったが、そもそもOuterとInnerの判別が主観によるもので、条件が曖昧すぎます。
という事で、自己交差している図形のテストを外して、テストを通過する事を目標に実装させました。

boost::geometry に足りないのは、make_valid に相当するルーチンだろうと思いまして、自社で作成したライブラリのアルゴリズムを make_valid として実装させる事にしました。
という感じで置いときます。 https://github.com/oki-hunes/boost_geometry_added

2026年3月9日月曜日

ArcGIS Experence Builder Deploy bug 忘備録

ArcGIS Experence Builder Deploymentには、
node -e "require('./server/src/middlewares/dev/apps/app-download.js').zipApp('0', 'app.zip', 'my_client_id');"
と書かれているが、出力された app.zip の config.json は、正しくないまま出力されている。

deploy の番号 '0' に対して
server/public/apps/0/resources/config/config.json
を、以下の場所にコピーしないと動作しない。
(cdn/3/ は、自動でインクリメントされる番号)
app/cdn/3/config.json
2026/3/10 ※コピー先を間違えていたので修正

2026年2月5日木曜日

EPSG忘備録

どうも、最近 JGD2011 以降の測地系を「日本測地系」と呼ぶらしい。そのせいで、Bessel1841の楕円体を地球の基準としていた日本の「旧日本測地系」のEPSG(SRID)が記載されなくなってきた。 しかしながら、仕事では、まだまだ旧日本測地系の成果を扱う事が多いので、メモっておかないと、わからなくなってしまう。

日本で取り扱う緯度経度座標には2種類あり

EPSG:4301: 旧日本測地系(ベッセル楕円体)
EPSG:4326: WGS84(GPSの標準)

となっている 旧日本測地系の平面直角座標
ESPG:30161: 旧日本測地系 平面直角座標1
ESPG:30162: 旧日本測地系 平面直角座標2
ESPG:30163: 旧日本測地系 平面直角座標3
ESPG:30164: 旧日本測地系 平面直角座標4
ESPG:30165: 旧日本測地系 平面直角座標5
ESPG:30166: 旧日本測地系 平面直角座標6
ESPG:30167: 旧日本測地系 平面直角座標7
ESPG:30168: 旧日本測地系 平面直角座標8
ESPG:30169: 旧日本測地系 平面直角座標9
ESPG:30170: 旧日本測地系 平面直角座標10
ESPG:30171: 旧日本測地系 平面直角座標11
ESPG:30172: 旧日本測地系 平面直角座標12
ESPG:30173: 旧日本測地系 平面直角座標13
ESPG:30174: 旧日本測地系 平面直角座標14
ESPG:30175: 旧日本測地系 平面直角座標15
ESPG:30176: 旧日本測地系 平面直角座標16

2021年4月20日火曜日

OpenCV warpPerspective を魔改造


ラスター変換で、図面の画像を縮小変換した場合に、Cubic や Lanczos4 といった周りの画素値から平均を求める手法だと、線が薄くなって消えてしまいます。縮小変換しても線を消さないようにするために、周りの画素値の最小値を取ってやれば良いと思い魔改造する事に・・・。

魔改造を施したのは OpenCV 4.5.2 のバージョンです。

/modules/imgproc/src/imgwarp.cpp
を改変していきます。尚、OpenCL や OpenVX など、コンパイル条件で多様なルーチンを選択できますが、最小限のコードのみを対象としました。
warpPerspectiveでは cv::INTER_AREA というフラグを指定しても cv::INTER_LINEAR というフラグに変換されて処理されていたので、cv::INTER_AREA が指定された場合に、線画用の縮小処理が実行されるという方針で実装しました

手っ取り早くBilinear のコードをベースに改造します。
係数配列を追加します。
static float AreasupTab_f[INTER_TAB_SIZE2][4][4];
static short AreasupTab_i[INTER_TAB_SIZE2][4][4];

初期化関数を追加します。不要なんですが、体裁だけ整えます。
static inline void interpolateAreasup( float x, float* coeffs )
{
    coeffs[0] = coeffs[1] = coeffs[2] = coeffs[3] = 0;
}

cv::INTER_AREAが指定された時に初期化されるコードを追加します
static const void* initInterTab2D( int method, bool fixpt )
{
    static bool inittab[INTER_MAX+1] = {false};
    float* tab = 0;
    short* itab = 0;
    int ksize = 0;
    if( method == INTER_LINEAR )
        tab = BilinearTab_f[0][0], itab = BilinearTab_i[0][0], ksize=2;
    else if( method == INTER_CUBIC )
        tab = BicubicTab_f[0][0], itab = BicubicTab_i[0][0], ksize=4;
    else if( method == INTER_LANCZOS4 )
        tab = Lanczos4Tab_f[0][0], itab = Lanczos4Tab_i[0][0], ksize=8;
        // --> ここから
    else if( method == INTER_AREA )
        tab = AreasupTab_f[0][0], itab = AreasupTab_i[0][0], ksize=4;
        // --> ここまで追加
    else
        CV_Error( CV_StsBadArg, "Unknown/unsupported interpolation type" );
    //...

係数配列を全部初期化する関数にもcv::INTER_AREA の初期化を追加します
static bool initAllInterTab2D()
{
    return  initInterTab2D( INTER_LINEAR, false ) &&
            initInterTab2D( INTER_LINEAR, true ) &&
            initInterTab2D( INTER_CUBIC, false ) &&
            initInterTab2D( INTER_CUBIC, true ) &&
            initInterTab2D( INTER_LANCZOS4, false ) &&
            initInterTab2D( INTER_LANCZOS4, true ) &&
            initInterTab2D( INTER_AREA, false ) && // この辺
            initInterTab2D( INTER_AREA, true );
}

処理の核心部を追加します
template<class CastOp, typename AT, int ONE>
static void remapAreasup( const Mat& _src, Mat& _dst, const Mat& _xy,
                          const Mat& _fxy, const void* _wtab,
                          int borderType, const Scalar& _borderValue )
{
    typedef typename CastOp::rtype T;
    typedef typename CastOp::type1 WT;
    Size ssize = _src.size(), dsize = _dst.size();
    const int cn = _src.channels();
    const AT* wtab = (const AT*)_wtab;
    const T* S0 = _src.ptr<T>();
    size_t sstep = _src.step/sizeof(S0[0]);
    T cval[CV_CN_MAX];
    CastOp castOp;

    for(int k = 0; k < cn; k++ )
        cval[k] = saturate_cast<T>(_borderValue[k & 3]);

    int borderType1 = borderType != BORDER_TRANSPARENT ? borderType : BORDER_REFLECT_101;

    unsigned width1 = std::max(ssize.width-3, 0), height1 = std::max(ssize.height-3, 0);

    if( _dst.isContinuous() && _xy.isContinuous() && _fxy.isContinuous() )
    {
        dsize.width *= dsize.height;
        dsize.height = 1;
    }

    for(int dy = 0; dy < dsize.height; dy++ )
    {
        T* D = _dst.ptr<T>(dy);
        const short* XY = _xy.ptr<short>(dy);
        const ushort* FXY = _fxy.ptr<ushort>(dy);

        for(int dx = 0; dx < dsize.width; dx++, D += cn )
        {
            int sx = XY[dx*2]-1, sy = XY[dx*2+1]-1;
            //const AT* w = wtab + FXY[dx]*16;
            if( (unsigned)sx < width1 && (unsigned)sy < height1 )
            {
                const T* S = S0 + sy*sstep + sx*cn;
                for(int k = 0; k < cn; k++ )
                {
                    WT sum = std::min<WT>(std::min<WT>(S[0],S[cn]), std::min<WT>(S[cn*2], S[cn*3]));
                    S += sstep;
                    sum = std::min<WT>(sum, std::min<WT>(std::min(S[0],S[cn]), std::min<WT>(S[cn*2], S[cn*3])));
                    S += sstep;
                    sum = std::min<WT>(sum, std::min<WT>(std::min(S[0],S[cn]), std::min<WT>(S[cn*2], S[cn*3])));
                    S += sstep;
                    sum = std::min<WT>(sum, std::min<WT>(std::min(S[0],S[cn]), std::min<WT>(S[cn*2], S[cn*3])));
                    S += 1 - sstep*3;
                    D[k] = static_cast<WT>(sum);
                }
            }
            else
            {
                int x[4], y[4];
                if( borderType == BORDER_TRANSPARENT &&
                    ((unsigned)(sx+1) >= (unsigned)ssize.width ||
                    (unsigned)(sy+1) >= (unsigned)ssize.height) ) {
                      continue;
                }

                if( borderType1 == BORDER_CONSTANT &&
                    (sx >= ssize.width || sx+4 <= 0 ||
                    sy >= ssize.height || sy+4 <= 0))
                {
                    for(int k = 0; k < cn; k++ )
                        D[k] = cval[k];
                    continue;
                }

                for(int i = 0; i < 4; i++ )
                {
                    x[i] = borderInterpolate(sx + i, ssize.width, borderType1)*cn;
                    y[i] = borderInterpolate(sy + i, ssize.height, borderType1);
                }

                for(int k = 0; k < cn; k++, S0++ )
                {
                  WT cv = cval[k], sum = std::numeric_limits<WT>::max();
                    for(int i = 0; i < 4; i++ )
                    {
                        int yi = y[i];
                        const T* S = S0 + yi*sstep;
                        if( yi < 0 )
                            continue;
                        if( x[0] >= 0 )
                            sum = std::min<WT>(sum, S[x[0]]);
                        if( x[1] >= 0 )
                            sum = std::min<WT>(sum, S[x[1]]);
                        if( x[2] >= 0 )
                            sum = std::min<WT>(sum, S[x[2]]);
                        if( x[3] >= 0 )
                            sum = std::min<WT>(sum, S[x[3]]);
                    }
                    D[k] = static_cast<WT>(sum);
                }
                S0 -= cn;
            }
        }
    }
}

remap に cv::INTER_AREA 時の処理を追加します。IPPとかは無視です。
void cv::remap( InputArray _src, OutputArray _dst,
                InputArray _map1, InputArray _map2,
                int interpolation, int borderType, const Scalar& borderValue )
{
    CV_INSTRUMENT_REGION();

    static RemapNNFunc nn_tab[] =
    {
        remapNearest<uchar>, remapNearest<schar>, remapNearest<ushort>, remapNearest<short>,
        remapNearest<int>, remapNearest<float>, remapNearest<double>, 0
    };

    static RemapFunc linear_tab[] =
    {
        remapBilinear<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, RemapVec_8u, short>, 0,
        remapBilinear<Cast<float, ushort>, RemapNoVec, float>,
        remapBilinear<Cast<float, short>, RemapNoVec, float>, 0,
        remapBilinear<Cast<float, float>, RemapNoVec, float>,
        remapBilinear<Cast<double, double>, RemapNoVec, float>, 0
    };

    static RemapFunc cubic_tab[] =
    {
        remapBicubic<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, short, INTER_REMAP_COEF_SCALE>, 0,
        remapBicubic<Cast<float, ushort>, float, 1>,
        remapBicubic<Cast<float, short>, float, 1>, 0,
        remapBicubic<Cast<float, float>, float, 1>,
        remapBicubic<Cast<double, double>, float, 1>, 0
    };

    static RemapFunc lanczos4_tab[] =
    {
        remapLanczos4<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, short, INTER_REMAP_COEF_SCALE>, 0,
        remapLanczos4<Cast<float, ushort>, float, 1>,
        remapLanczos4<Cast<float, short>, float, 1>, 0,
        remapLanczos4<Cast<float, float>, float, 1>,
        remapLanczos4<Cast<double, double>, float, 1>, 0
    };

    // --> ここから
    static RemapFunc areasup_tab[] = 
    {
        remapAreasup<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, short, INTER_REMAP_COEF_SCALE>, 0,
        remapAreasup<Cast<float, ushort>, float, 1>,
        remapAreasup<Cast<float, short>, float, 1>, 0,
        remapAreasup<Cast<float, float>, float, 1>,
        remapAreasup<Cast<double, double>, float, 1>, 0
    };
    // --> ここまで追加
  
    CV_Assert( !_map1.empty() );
    CV_Assert( _map2.empty() || (_map2.size() == _map1.size()));

    CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
               ocl_remap(_src, _dst, _map1, _map2, interpolation, borderType, borderValue))

    Mat src = _src.getMat(), map1 = _map1.getMat(), map2 = _map2.getMat();
    _dst.create( map1.size(), src.type() );
    Mat dst = _dst.getMat();


    CV_OVX_RUN(
        src.type() == CV_8UC1 && dst.type() == CV_8UC1 &&
        !ovx::skipSmallImages<VX_KERNEL_REMAP>(src.cols, src.rows) &&
        (borderType& ~BORDER_ISOLATED) == BORDER_CONSTANT &&
        ((map1.type() == CV_32FC2 && map2.empty() && map1.size == dst.size) ||
         (map1.type() == CV_32FC1 && map2.type() == CV_32FC1 && map1.size == dst.size && map2.size == dst.size) ||
         (map1.empty() && map2.type() == CV_32FC2 && map2.size == dst.size)) &&
        ((borderType & BORDER_ISOLATED) != 0 || !src.isSubmatrix()),
        openvx_remap(src, dst, map1, map2, interpolation, borderValue));

    CV_Assert( dst.cols < SHRT_MAX && dst.rows < SHRT_MAX && src.cols < SHRT_MAX && src.rows < SHRT_MAX );

    if( dst.data == src.data )
        src = src.clone();

    // --> 使われていない INTER_AREA を利用するので、以下コメントアウトして有効化
    //if( interpolation == INTER_AREA )
    //    interpolation = INTER_LINEAR;

    int type = src.type(), depth = CV_MAT_DEPTH(type);

#if defined HAVE_IPP && !IPP_DISABLE_REMAP
    CV_IPP_CHECK()
    {
        if ((interpolation == INTER_LINEAR || interpolation == INTER_CUBIC || interpolation == INTER_NEAREST) &&
                map1.type() == CV_32FC1 && map2.type() == CV_32FC1 &&
                (borderType == BORDER_CONSTANT || borderType == BORDER_TRANSPARENT))
        {
            int ippInterpolation =
                interpolation == INTER_NEAREST ? IPPI_INTER_NN :
                interpolation == INTER_LINEAR ? IPPI_INTER_LINEAR : IPPI_INTER_CUBIC;

            ippiRemap ippFunc =
                type == CV_8UC1 ? (ippiRemap)ippiRemap_8u_C1R :
                type == CV_8UC3 ? (ippiRemap)ippiRemap_8u_C3R :
                type == CV_8UC4 ? (ippiRemap)ippiRemap_8u_C4R :
                type == CV_16UC1 ? (ippiRemap)ippiRemap_16u_C1R :
                type == CV_16UC3 ? (ippiRemap)ippiRemap_16u_C3R :
                type == CV_16UC4 ? (ippiRemap)ippiRemap_16u_C4R :
                type == CV_32FC1 ? (ippiRemap)ippiRemap_32f_C1R :
                type == CV_32FC3 ? (ippiRemap)ippiRemap_32f_C3R :
                type == CV_32FC4 ? (ippiRemap)ippiRemap_32f_C4R : 0;

            if (ippFunc)
            {
                bool ok;
                IPPRemapInvoker invoker(src, dst, map1, map2, ippFunc, ippInterpolation,
                                        borderType, borderValue, &ok);
                Range range(0, dst.rows);
                parallel_for_(range, invoker, dst.total() / (double)(1 << 16));

                if (ok)
                {
                    CV_IMPL_ADD(CV_IMPL_IPP|CV_IMPL_MT);
                    return;
                }
                setIppErrorStatus();
            }
        }
    }
#endif

    RemapNNFunc nnfunc = 0;
    RemapFunc ifunc = 0;
    const void* ctab = 0;
    bool fixpt = depth == CV_8U;
    bool planar_input = false;

    if( interpolation == INTER_NEAREST )
    {
        nnfunc = nn_tab[depth];
        CV_Assert( nnfunc != 0 );
    }
    else
    {
        if( interpolation == INTER_LINEAR )
            ifunc = linear_tab[depth];
        else if( interpolation == INTER_CUBIC ){
            ifunc = cubic_tab[depth];
            CV_Assert( _src.channels() <= 4 );
        }
        else if( interpolation == INTER_LANCZOS4 ){
            ifunc = lanczos4_tab[depth];
            CV_Assert( _src.channels() <= 4 );
        }
        // --> ここから
        else if( interpolation == INTER_AREA ) {
            ifunc = areasup_tab[depth];
            CV_Assert( _src.channels() <= 4 );
        }
        // --> ここまで追加
        else
            CV_Error( CV_StsBadArg, "Unknown interpolation method" );
        CV_Assert( ifunc != 0 );
        ctab = initInterTab2D( interpolation, fixpt );
    }

    const Mat *m1 = &map1, *m2 = &map2;

    if( (map1.type() == CV_16SC2 && (map2.type() == CV_16UC1 || map2.type() == CV_16SC1 || map2.empty())) ||
        (map2.type() == CV_16SC2 && (map1.type() == CV_16UC1 || map1.type() == CV_16SC1 || map1.empty())) )
    {
        if( map1.type() != CV_16SC2 )
            std::swap(m1, m2);
    }
    else
    {
        CV_Assert( ((map1.type() == CV_32FC2 || map1.type() == CV_16SC2) && map2.empty()) ||
            (map1.type() == CV_32FC1 && map2.type() == CV_32FC1) );
        planar_input = map1.channels() == 1;
    }

    RemapInvoker invoker(src, dst, m1, m2,
                         borderType, borderValue, planar_input, nnfunc, ifunc,
                         ctab);
    parallel_for_(Range(0, dst.rows), invoker, dst.total()/(double)(1<<16));
}

成果
 変換元画像

 cv::INTER_LINEAR で変換

 魔改造 cv::INTER_AREA で変換

2021年2月3日水曜日

Proj4 C++ 座標変換忘備録

Proj4 ver 6.3.1 で座標変換するコードの習作
#include <proj.h>
#include <iostream>

int main() {

  PJ_CONTEXT* C;
  PJ*         P;
  PJ*         P_for_GIS;

  C = proj_context_create();
  P = proj_create_crs_to_crs( C, "EPSG:3857", "EPSG:6679", NULL );
  if( !P ) {
    std::cerr << "Fail to create proj crs conversion" << std::endl;
    return 1;
  }
  P_for_GIS = proj_normalize_for_visualization(C,P);
  if( !P_for_GIS ) {
    std::cerr << "Fail to normalize proj" << std::endl;
    return 1;
  }
  proj_destroy(P);
  P = P_for_GIS;
  
  PJ_COORD sp = proj_coord(15734773.4,5323579.9,0,0);
  
  PJ_COORD dp = proj_trans(P, PJ_FWD, sp);
  std::cout << std::fixed << sp.v[0] << "," << sp.v[1] << std::endl;
  std::cout << " to ";
  std::cout << dp.v[0] << "," << dp.v[1] << std::endl;
  
  proj_destroy(P);
  proj_context_destroy(C);
  
  return 0;
}
参考: https://github.com/OSGeo/PROJ/blob/master/examples/pj_obs_api_mini_demo.c

2016年12月15日木曜日

tile splash

はじめに

ウェブ上で地図の配信を行う手法は、タイルに分割された地図画像をajaxにより書き換える方法と、ベクトルデータを全部・あるいは部分的に読み込む手法が主でした。近年、ベクトルデータをタイル分割し、ajax により書き換える方法が登場しました。それが VectorTileです。
 このVectorTileの配信を行うオープンソースのモジュールのうちのひとつがTileSplashです。TileSplash は node.js という javascript によるWebサーバのエンジン用ライブラリとして書かれています。プロジェクトは、GitHubで下記にて公開されています。

 https://github.com/faradayio/tilesplash

実行に必要なもの

TileSplash を実行させるには、以下のモジュールが必要です。 node.js npm PostgreSQL server PostGIS インストール方法については、割愛します。

TileSplash 環境の構築

コマンド・プロンプトにて、npm を利用して SplashTile をインストールします。
$ npm install tilesplash
尚、-g オプションを使用すると、グローバル環境にモジュールがインストールされますが、ここでは、コマンドを実行したディレクトリ上にモジュールがインストールされるようにします。
 PostgreSQL Server に図形テーブルを構築します。注意点として、図形の座標系は、EPSG4326 でないと動作しません。また、ST_Transform(geom,4326) と言った書き方をしてもエラーになるので注意してください。既存のgeometryフィールドが異なるSRIDで作成している場合は、以下のようにフィールドを追加すると良いでしょう。
 alter table target_table add column geom2 geometry(‘MultiPolygon’, 4326);
    update table target_table set geom2 = ST_Transform(geom,4326); 
エディタで、sptile.js を作成し編集します。
var Tilesplash = require('tilesplash');
// username に postgresql へ接続するユーザ名を指定
// localhost にpostgresql server ホスト名を指定
// dbname に postgresql のデータベース名を指定
var app = new Tilesplash('postgres://username@localhost/dbname);

// layer_name はレイヤ名で、ajax によるリクエストの一部として扱われる
// table_name は図形テーブル名
// geom は図形フィールド名
// ここではシンプリファイをかけていますが、tile パラメータの x 値により
// シンプリファイをかけるパラメータを変更する必要があると思われます
app.layer('layer_name', function(tile, render){
  render('SELECT ST_AsGeoJSON(ST_Simplify(geom,0.005)) as the_geom_geojson’
     ‘FROM target_table WHERE ST_Intersects(St_Simplify(geom,0.005), !bbox_4326!)');
});
// ポート3000にて、Webサーバを起動します
app.server.listen(3000);
このまま稼働させると、レスポンスが javascript になるため、違うサイト間でのリクエストで CORS(Cross Origin Resource Sharing)の制約にひっかかって動作しない可能性があります。具体的には「No 'Access-Control-Allow-Origin' header is present」といったエラーがブラウザにより返されます。これを解消するには、以下のモジュールを利用します。
http://www.slideshare.net/kitfactory/web-api-34814937
$npm install corser

でモジュールをインストールし、
node_mojules/tilesplash/lib/index.js
ファイルに以下の修正を施します。
  this.server.use(pgMiddleware(dbOptions));

  var corser = require("corser");      // 以下の2行を CORS 対策として挿入
  this.server.use(corser.create());

  this.cacheOptions = cacheOptions || {};
  this._cache = new Caching(cacheType || 'memory', cacheOptions);
これは、レスポンス・ヘッダに
  Access-Control-Allow-Origin: *

を付加する事になります。

クライアント実装例

ローカルホストの node サーバに対して、アクセスする例です。
<!DOCTYPE html>
<html>
<head>
<meta charset="UTF-8">
<title>GSI Tiles on OpenLayers 3</title>
<link rel="stylesheet" href="http://openlayers.org/en/v3.10.1/css/ol.css" type="text/css">
<link rel="stylesheet" href="http://maxcdn.bootstrapcdn.com/bootstrap/3.2.0/css/bootstrap.min.css">
<script src="http://ajax.googleapis.com/ajax/libs/jquery/1.11.1/jquery.min.js"></script>
<script src="http://maxcdn.bootstrapcdn.com/bootstrap/3.2.0/js/bootstrap.min.js"></script>
<script src="http://openlayers.org/en/v3.10.1/build/ol.js" type="text/javascript"></script>


<style>
  body {padding: 0; margin: 0}
  html, body, #map {height: 100%; width: 100%;}
</style>
</head>
<body>
  <div id="page">
    <div id="head"></div>
    <div id="main">  
      <div id="map" style="float:right"; width:640px; margin:0; padding:0; ></div>
      <div id="left" style="float:left; width:300px; margin:0; padding:0;"></div>
    </div>
  </div>       
<script>
var map = new ol.Map({
  target: "map",
  renderer: ['canvas', 'dom'],
  layers: [
    // 地理院タイル・レイヤ
    new ol.layer.Tile({
      source: new ol.source.XYZ({
        attributions: [
          new ol.Attribution({
            html: "<a href='http://maps.gsi.go.jp/development/ichiran.html' target='_blank'>地理院タイル</a>"
          })
        ],
        url: "http://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png",
        projection: "EPSG:3857"
      })
    }),
    
    // topoJson レイヤ (splashtile から引く)
    new ol.layer.Vector({
      source: new ol.source.TileVector({
        format: new ol.format.TopoJSON(),
        projection: 'EPSG:3857',
        tileGrid: new ol.tilegrid.createXYZ({
          maxZoom: 19
        }),
        url: 'http://localhost:3000/' +
          'test_layer/{z}/{x}/{y}.topojson'
        }),
       style: new ol.style.Style({
         fill: new ol.style.Fill({ color: '#9db9e8' }),
         stroke: new ol.style.Stroke({ color: '#FF0000' })
       })
    })
  ],
  controls: ol.control.defaults({
    attributionOptions: ({
      collapsible: false
    })
  }),
  view: new ol.View({
    projection: "EPSG:3857",
    center: ol.proj.transform([138.7313889, 35.3622222], "EPSG:4326", "EPSG:3857"),
    maxZoom: 18,
    zoom: 5
    })
});
</script>
</body>
</html>
実行すると、こんな感じになります。 file:/// 上から実行すると、CORSにひっかかります。chrome ではなく、safari で実行しました。 CORS対策の部分は、もうちょい真面目にやらないとダメかもしれません。

2014年11月6日木曜日

GDAL-1.11.1 における ShapeDriver の日本語設定

前回「GDAL-1.10.1 の日本語環境ではまった」からの続きです。 うまく行ったり、行かなかったりは、文字コードが UTF-8 だったり SJIS だったりと、判別のしようが無い Shape が混在していたのが大きな原因でした。 それで、定説の 環境変数
set SHAPE_ENCODING=LDID/19 
ですが、そこを解釈するコードは、このようになっております。
CPLString OGRShapeLayer::ConvertCodePage( const char *pszCodePage )

{
    CPLString osEncoding;
  osEncoding.Clear();

    if( pszCodePage == NULL )
        return osEncoding;

    if( EQUALN(pszCodePage,"LDID/",5) )
    {
        int nCP = -1; // windows code page. 

        //http://www.autopark.ru/ASBProgrammerGuide/DBFSTRUC.HTM
        switch( atoi(pszCodePage+5) )
        {
          case 1: nCP = 437;      break;
          case 2: nCP = 850;      break;
          case 3: nCP = 1252;     break;
          case 4: nCP = 10000;    break;
          case 8: nCP = 865;      break;
          case 10: nCP = 850;     break;
          case 11: nCP = 437;     break;
          case 13: nCP = 437;     break;
          case 14: nCP = 850;     break;
          case 15: nCP = 437;     break;
          case 16: nCP = 850;     break;
          case 17: nCP = 437;     break;
          case 18: nCP = 850;     break;
          case 19: nCP = 932;     break;
          case 20: nCP = 850;     break;
          case 21: nCP = 437;     break;
          case 22: nCP = 850;     break;
          case 23: nCP = 865;     break;
          case 24: nCP = 437;     break;
          case 25: nCP = 437;     break;
          case 26: nCP = 850;     break;
          case 27: nCP = 437;     break;
          case 28: nCP = 863;     break;
          case 29: nCP = 850;     break;
          case 31: nCP = 852;     break;
          case 34: nCP = 852;     break;
          case 35: nCP = 852;     break;
          case 36: nCP = 860;     break;
          case 37: nCP = 850;     break;
          case 38: nCP = 866;     break;
          case 55: nCP = 850;     break;
          case 64: nCP = 852;     break;
          case 77: nCP = 936;     break;
          case 78: nCP = 949;     break;
          case 79: nCP = 950;     break;
          case 80: nCP = 874;     break;
          case 87: return CPL_ENC_ISO8859_1;
          case 88: nCP = 1252;     break;
          case 89: nCP = 1252;     break;
          case 100: nCP = 852;     break;
          case 101: nCP = 866;     break;
          case 102: nCP = 865;     break;
          case 103: nCP = 861;     break;
          case 104: nCP = 895;     break;
          case 105: nCP = 620;     break;
          case 106: nCP = 737;     break;
          case 107: nCP = 857;     break;
          case 108: nCP = 863;     break;
          case 120: nCP = 950;     break;
          case 121: nCP = 949;     break;
          case 122: nCP = 936;     break;
          case 123: nCP = 932;     break;
          case 124: nCP = 874;     break;
          case 134: nCP = 737;     break;
          case 135: nCP = 852;     break;
          case 136: nCP = 857;     break;
          case 150: nCP = 10007;   break;
          case 151: nCP = 10029;   break;
          case 200: nCP = 1250;    break;
          case 201: nCP = 1251;    break;
          case 202: nCP = 1254;    break;
          case 203: nCP = 1253;    break;
          case 204: nCP = 1257;    break;
          default: break;
        }

        if( nCP != -1 )
        {
            osEncoding.Printf( "CP%d", nCP );
            return osEncoding;
        }
    }

    // From the CPG file
    // http://resources.arcgis.com/fr/content/kbase?fa=articleShow&d=21106
    
    if( (atoi(pszCodePage) >= 437 && atoi(pszCodePage) <= 950)
        || (atoi(pszCodePage) >= 1250 && atoi(pszCodePage) <= 1258) )
    {
        osEncoding.Printf( "CP%d", atoi(pszCodePage) );
        return osEncoding;
    }
    if( EQUALN(pszCodePage,"8859",4) )
    {
        if( pszCodePage[4] == '-' )
            osEncoding.Printf( "ISO-8859-%s", pszCodePage + 5 );
        else
            osEncoding.Printf( "ISO-8859-%s", pszCodePage + 4 );
        return osEncoding;
    }
    if( EQUALN(pszCodePage,"UTF-8",5) )
        return CPL_ENC_UTF8;

    // try just using the CPG value directly.  Works for stuff like Big5.
    return pszCodePage;
}
これのディスカッションは、ESRI社のフォーラムで話題になっていました。
コード中の解説ページhttp://www.autopark.ru/ASBProgrammerGuide/DBFSTRUC.HTMが詳しいです。
CP932 に相当するケース文の値は、0x7B でなければなりませんから、10進数にしますと、123になります。上記コードと照らしあわせても整合がとれていますよね? という事で、LDID/19 でもOKなのですが、
set SHAPE_ENCODING=LDID/123 
を設定すると、CP932 に変換してくれる。
もしくは、最初から CP932 を指定していれば、LDID/の変換をパスして、上記関数を通った後の値が評価されるので
set SHAPE_ENCODING=CP932 
でもOKだったという事になるのでしょうか? 拡張子 DBF ファイルの 29バイト目に、case 文の値か、UTF-8なら -1 == 65535 == 0xFFFF を設定するように改善していけば良さそうですよね?
と、思ったら、-1 は CPACP 扱い? int は 32bit 以上で、65535 == 0xFFFF の解釈で考えればOKでしょうか? 同日追記:  実は SHAPE_ENCODING= を設定しない方が、うまく動く場合もあるようです。どうも、前回DBFファイルに設定があった場合は、SHAPE_ENCODINGは使用されないような事を書きましたが、その記述が間違いだったみたいです。そして、SHAPE_ENCODING=UTF-8 と設定してみましたが、挙動が違うような・・・。もうちょい、突っ込んで調べてみます・・・

2013年4月3日水曜日

GeoTIFFにやられまくり

TIFFTAG_GDAL_NODATA でハマった事を書きましたが、続きがあります。
GeoTIFF では、GDAL_NODATA に設定された値がアルファ値0、つまり透過として扱われる事になります。
RGB3バンドあった場合に、R == 255 && G == 255 && B==255 なら透過扱いなんかなぁ?と思ってたんですが、どうも実装する人によって解釈が違ってて、RGB(255,0,0) が透過色に設定されてたりしました。こうなってくると、RGB(254,255,255) は透過扱いなのか、透過扱いでないのか?これすら怪しいです。
じゃあ、インデックス・カラーの場合は、どうか?というと、カラーインデックス番号 255 が透過扱いでした。いや、もうカオスで、嫌になってきます。
GDAL_NODATA が設定されてるんだから、もちろん、タイル内のブロック全てに、この値が適応されてるんですよね?と思ったら、画像高を越えた領域には、設定されていませんでした。お陰で、黒い領域が出現して何事かと思いましたよ。

関係ないですけど、なんか、ブログの RSS フィードがおかしくなって出力されないですし、Blogger サービスもいつまであるかわからないので、引越しを検討した方がいいのかなぁ?なんて考えてます。

2012年7月12日木曜日

QGis Python Plugin と日本語の周辺1

QGIS の python plugin を書いていて、やはり日本語の処理に四苦八苦したので、判明した事をチビチビと書いていこうと思います。ノウハウが溜まる度に書き足すつもりなので、その1としました。どれぐらいの頻度で続くのかは、書いてる本人にもわかりません。

 文字コードの国際化に関する状況は、刻々と変化する可能性が高いので、環境を明記しときます。
  • Windows 7 (x64)
  • GDAL 1.9.1
  • Quatom GIS Lisboa 1.8.0

 Windows 環境において、GDAL/OGR あたりで、日本語のファイル名を使おうと思ったら、環境変数「GDAL_FILENAME_IS_UTF8」に「NO」と設定しないと、ogr2ogr 等のコマンドでファイル名が不正だとエラーにされてしまう現状があります。

 それなら、思い切ってコントロールパネルの「システム環境設定」に、GDAL_FILENAME_IS_UTF8 = NO と設定してしまえば良いじゃないか?と設定したら、ものの見事に嵌りました。なんと、QGISでshapefileを開こうとするとエラーになってしまいます。なんてこったい!QGIS内部では、UTF8でファイル名を扱っている雰囲気が漂ってます。これに関しては、現段階では確認できていません。

 プラグインのコードで、QFileDialog からファイル名を取得するコードを見てみましょう

    fileName = QFileDialog.getSaveFileName(None,QString.fromLocal8Bit("Select a file:"), "", "*.tapgis")
    if fileName.isNull():
      return
    fname = unicode(fileName.toLocal8Bit(),'mbcs')
  はい、ウィンドウズでは、ファイル名はOEMの文字コードセットになりますので、multi byte char set  つまり、SJIS(CP932)になっております。python なんて、ほとんど知らねぇよだったので、unicode という名前の関数が、文字コード・セットを指定して文字列を初期化する関数だと判るまでに丸一日費やしました。
  さて、こいつを ogr2ogr のコマンドに引き渡すとします。

  実はファイル名に半角スペースが混じっている場合には、”path" というようにダブルクォーテーションでパスを括ってやらなければなりません。逆に半角スペースが混じっていない場合には、path というようにダブルクォーテーションで囲ってはいけません(何この仕様・・・くたばれ)。という事情は、すっ飛ばします。
  # ogrvars には "-overwrite" 等のオプション文字列が渡される
  def CallRWTools(self, infile, outfile, ogrvars):
    import os
    try:
      import subprocess
      path = str(QgsApplication.prefixPath()).rsplit('/', 2)[0]
      pcmd = path + '\\bin\\ogr2ogr';
      fmt = '%s -f "tapgis" %s %s %s'
      cmd = fmt % (pcmd, outfile, infile, ogrvars)
      # ファイル名がUTF8ではないという状態に変更
      os.environ["GDAL_FILENAME_IS_UTF8"] = "NO"
      subprocess.call(cmd.encode('mbcs'), shell=True)
      #cmd3 = "echo " + cmd + " > " + outfile + ".txt"
      #subprocess.call(cmd3.encode('mbcs'), shell=True)
      #cmd2 = "notepad " + outfile + ".txt"
      #subprocess.call(cmd2.encode('mbcs'), shell=True)
    finally:
      # ファイル名がUTF8であるという状態に戻す
      os.environ["GDAL_FILENAME_IS_UTF8"] = "YES"
  はい、QGIS 内部では、GDAL_FILENAME_IS_UTF8 = YES を前提としているので、コマンドを呼び出す間だけ、GDAL_FILENAME_IS_UTF8 = NO に変えてコールします。あと、WIndows シェルでは mbcs なので、encode('mbcs') してやらないと、正しいコマンド引数になりません。

  ここまで書いても、日本語の問題は残ります。今度は、ogr2ogr のオプション文字列が mbcs のままドライバに渡されるためです。自分の書いたドライバでは、sqlite3_open という関数をコールしているので、OGRDriver の OGRDataSource::Open, OGRDataSource::Create という関数内で渡されたファイル名が mbcs なので utf8 へ文字コード変換してやる必要がありました。   また、OGRDataSource::CreateLayer に渡されるレイヤ名は、変換元のデータソースのOGRDriver の実装に依存している状態なので、どんな文字コードが渡されるのか、未知数です。shape driver の実装では、-nln オプションでレイヤ名を指定しない場合には、ファイル名からレイヤ名を生成しているので、mbcs の文字コードが来ました。注意しなければならないのは、これは Windows の日本語OS 独特の状況においてのみ発生します。Windows の Laten1なOSの環境では、おそらく Laten1の文字コードが来るという結果になるでしょう。-nln オプションも生の文字コードが渡されるので、同様の結果になります。

  このように、locale でなんとかしようという発想は、破綻するので utf8 で統一的に扱う必要性が高いと思います。いちいちロケール指定してたら大変です。ウィンドウズの場合は、CPACP を有効活用する方向が望ましいと思います。 2012/07/19 追記: その2に続く 2012/09/05 追記:結局、ウィンドウズ環境では、GDAL_FILENAME_IS_UTF8 = YES を前提にコードが書かれてあるので、その1のアプローチはダメダメでした。

2012年7月2日月曜日

FOSS4G HOKKAIDO 2012 で発表しました

FOSS4G HOKKAIDO 2012 で発表してきました。資料は、こちらになります。  15分間に3つもテーマを入れるのは、詰め込み過ぎたかもしれません。これでも、まだネタは有って、テーマは絞りました。デモも交えようかと考えていたのですが、さすがに時間的にちょっと無理で、断念しました。  最初のDynamicTileというやつですが、会社のデモページの機能選択で「地図色塗り解析」というやつで使用されています。実際には、php側で画像ファイルはディスク上にキャッシュをかけていて、1日に1回ぐらいのリセットをかけてます。他には、機能選択の「レイヤ設定」で画面左に出てくる「表示レイヤ選択」というやつも、該当します。こいつも、レイヤの組み合わせで合成したやつをディスク上にキャッシュをかける仕組みになっています。  今年の春に参加した DevFestX の感想でも、ぼやいているように、FusionTable というやつと似てます。時期的には、開発されたのは同時期ぐらいだと思われます。  2番目のタブレット型GISは、現在進行形で開発しているものです。道内の自治体に1カ所だけ納品しています。一般売りできるように、改良を重ねてます。ビデオを制作してもらったので、自分が解説するより、わかりやすいかな?と思って、そのまま流しました。宣伝になってしまいますが、3番目の発表に続く文脈なので、ご容赦ください(^^;  3番目のGDAL/OGRの話は、先ほど発表された基盤地図対応GDAL/OGRを見て、ここら辺の話をできる人間って、そんなにいないのかな?と思って、まとめてみました。若干会場を置き去りにしてしまったかもしれません...。  発表中の   Android 用 Proj4j は、こちら   GDAL/OGR スタートアップ・キットは、こちら  になります。GDAL/OGR の方は、QGis を C:¥OSGEO4W という場所にインストールした事を前提にしていて、検索パスを書いてます。パスが違ってうまく行かない場合には、cmake_modules/FindOGR.cmake というテキストファイルを編集して、試してみてください。CMakeLists.txt の冒頭に cmake のコマンドの使い方が書いてあります。ドキュメント無しは、ちょっと味気なかったですね...。  何かあったら、Google+ 沖 観行 まで、メンション飛ばしてください。あと、OSGeoJP の ML にも登録するようにしますので、よろしくお願いします。  最後に、北海道地図の朝日さま、原田さま、デジタル北海道の三好さま、日頃仕事でお世話になっている北海道地図さま、OSGeoJP 関係者の皆様、ありがとうございました。