Try to make a CS 3D tile from the Geographical Survey tile
Overview
The existence of a map "CS 3D map" with emphasized micro-topography, which was recently devised by Nagano Forestry Research Center in connection with work. I found out. The CS 3D map is created by synthesizing the slope map and the curvature map (Curvature), and is a 3D map that makes it easy to read the terrain such as unevenness. Various tools such as QGIs plugin are already provided, but for studying, I tried to create CS 3D tile from Geographical Survey tile. It was.
How to make
Since the tilt map required to create a CS 3D map is provided by the Geographical Survey Institute tile, I thought that if I create a curvature map from the DEM tile, I can create a CS 3D map. Therefore, the creation procedure is as follows.
- Obtain the slope map tile from the Geographical Survey Institute tile.
- Obtain DEM tile from Geographical Survey tile
- Calculate the curvature from DEM and create a curvature diagram tile
- Multiply the pixels of the tilt map tile and the curvature map tile to generate a CS solid map tile
Source code
The source code is uploaded on Github. -> https://github.com/termat/CSMapTile Please see there for details. The CS 3D tile is generated by the TileDB class, and the main processing contents are as follows.
1. CS 3D tile generation
TileDB.java (excerpt)
private static String slope="https://cyberjapandata.gsi.go.jp/xyz/slopemap";
private static String dem14="https://cyberjapandata.gsi.go.jp/xyz/dem_png";
private static String dem15a="https://cyberjapandata.gsi.go.jp/xyz/dem5a_png";
private static String dem15b="https://cyberjapandata.gsi.go.jp/xyz/dem5b_png";
private static final int P8=(int)Math.pow(2,8);
private static final int P16=(int)Math.pow(2,16);
private static final int P23=(int)Math.pow(2,23);
private static final int P24=(int)Math.pow(2,24);
private static final double U=0.01;
private static double GEO_R=6378137;
private static int[] col=new int[]{255,232,197};
/*
*CS tile image acquisition
*/
public BufferedImage getCSImage(int zoom,int tileX,int tileY)throws IOException{
String param="/"+zoom+"/"+tileX+"/"+tileY;
BufferedImage s_img=ImageIO.read(new URL(slope+param+".png "));
s_img=getColorTransSlope(s_img);
BufferedImage d_img=createDemTile(zoom,tileX,tileY);
BufferedImage c_img=getCurve(d_img,zoom);
BufferedImage img=mul(c_img,s_img);
return img;
}
/*
*Convert the tilt map to orange
*/
private static BufferedImage getColorTransSlope(BufferedImage img){
int w=img.getWidth();
int h=img.getHeight();
BufferedImage ret=new BufferedImage(w,h,img.getType());
for(int i=0;i<w;i++){
for(int j=0;j<h;j++){
float cv=(float)(img.getRGB(i, j)&0x0000ff)/255f;
ret.setRGB(i, j, new Color((int)(col[0]*cv),(int)(col[1]*cv),(int)(col[2]*cv)).getRGB());
}
}
return ret;
}
/*
*Generate 258 x 258 size Dem tiles
*/
private static BufferedImage createDemTile(int zoom,int tileX,int tileY)throws IOException{
BufferedImage bc=getDemImage(zoom,tileX,tileY);
BufferedImage ret=new BufferedImage(bc.getWidth()+2,bc.getHeight()+2,bc.getType());
Graphics2D g=ret.createGraphics();
g.drawImage(bc, 1, 1, comp);
try{
BufferedImage bu=getDemImage(zoom,tileX,tileY-1);
g.drawImage(bu, 1, -255, comp);
}catch(IOException e){}
try{
BufferedImage bd=getDemImage(zoom,tileX,tileY+1);
g.drawImage(bd, 1, 257, comp);
}catch(IOException e){}
try{
BufferedImage bl=getDemImage(zoom,tileX-1,tileY);
g.drawImage(bl, -255, 1, comp);
}catch(IOException e){}
try{
BufferedImage br=getDemImage(zoom,tileX+1,tileY);
g.drawImage(br, 257, 1, comp);
}catch(IOException e){}
g.dispose();
return ret;
}
/*
*Obtained Dem tile from Geographical Survey Institute
*/
private static BufferedImage getDemImage(int zoom,int tileX,int tileY)throws IOException{
String param="/"+zoom+"/"+tileX+"/"+tileY;
BufferedImage d_img=null;
if(zoom<=14){
d_img=ImageIO.read(new URL(dem14+param+".png "));
}else{
try{
d_img=ImageIO.read(new URL(dem15a+param+".png "));
}catch(IOException e){
d_img=ImageIO.read(new URL(dem15b+param+".png "));
}
}
return d_img;
}
/*
*Image multiplication
*/
private static BufferedImage mul(BufferedImage im1,BufferedImage im2){
int w=im1.getWidth();
int h=im1.getHeight();
BufferedImage ret=new BufferedImage(w,h,BufferedImage.TYPE_INT_RGB);
for(int i=0;i<w;i++){
for(int j=0;j<h;j++){
int rgb1=im1.getRGB(i, j);
int rgb2=im2.getRGB(i, j);
ret.setRGB(i, j, mulRGB(rgb1,rgb2));
}
}
return ret;
}
/*
*RGB multiplication
*/
private static int mulRGB(int rgb1,int rgb2){
float r1=(float)(rgb1 >> 16 & 0xff)/255f;
float g1=(float)(rgb1&0x00ff00 >> 8 & 0xff)/255f;
float b1=(float)(rgb1&0x0000ff & 0xff)/255f;
float r2=(float)(rgb2 >> 16 & 0xff)/255f;
float g2=(float)(rgb2&0x00ff00 >> 8 & 0xff)/255f;
float b2=(float)(rgb2&0x0000ff & 0xff)/255f;
return new Color(r1*r2,g1*g2,b1*b2).getRGB();
}
/*
*Curvature image acquisition
*/
public static BufferedImage getCurve(BufferedImage dem,int zoom){
BufferedImage ret=new BufferedImage(dem.getWidth(),dem.getHeight(),BufferedImage.TYPE_INT_RGB);
double[][] dd=new double[dem.getWidth()][dem.getHeight()];
for(int i=0;i<dd.length;i++){
for(int j=0;j<dd[i].length;j++){
dd[i][j]=getZ(dem.getRGB(i, j));
}
}
double ll=getL(zoom);
for(int i=1;i<dd.length-1;i++){
for(int j=1;j<dd[i].length-1;j++){
double[][] p=new double[][]{
{dd[i-1][j-1],dd[i-1][j],dd[i-1][j+1]},
{dd[i][j-1],dd[i][j],dd[i][j+1]},
{dd[i+1][j-1],dd[i+1][j],dd[i+1][j+1]}};
double cu=getCurveVal(p,ll);
Color col=grad1.getColor(range.getNormalValue(cu));
ret.setRGB(i, j, col.getRGB());
}
}
return ret.getSubimage(1, 1, 256, 256);
}
/*
*Elevation acquisition
*/
private static double getZ(int intColor){
Color c=new Color(intColor);
int r=c.getRed();
int g=c.getGreen();
int b=c.getBlue();
int x=r*P16+g*P8+b;
if(x<P23){
return U*(double)x;
}else if(x>P23){
return 0;
}else{
return U*(double)(x-P24);
}
}
/*
*Zoom level 1 pixel length(m)
*/
private static double getL(int zoom){
return 2*GEO_R*Math.PI/256/Math.pow(2, zoom);
}
/*
*Curvature calculation
*/
private static double getCurveVal(double[][] p,double ll){
double ex1=p[1][0]+p[1][2]-2*p[1][1];
double ex2=p[0][1]+p[2][1]-2*p[1][1];
return (ex1+ex2)/(ll*ll)*100;
}
■ Geospatial Information Authority of Japan Inclination tile
■ Geospatial Information Authority of Japan DEM tile
■ Curvature diagram tile calculated from DEM tile
■ Generated CS 3D tile
2. Landslide terrain tile generation
Since it is a great deal, landslide topography information from Landslide Topography Distribution Map WMS Service published by the National Research Institute for Earth Science and Disaster Prevention I decided to acquire and tile it and overlay it with the CM 3D drawing tile. Acquisition and tiling of landslide topography information is handled by the JShisWMS class.
JShisWMS.java
import java.awt.image.BufferedImage;
import java.net.URL;
import javax.imageio.ImageIO;
public class JShisWMS {
private static final String baseURL="http://www.j-shis.bosai.go.jp/map/wms/landslide?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&";
private static final double L=(180.0/Math.PI)*Math.asin(Math.tanh(Math.PI));
/*
* http://www.j-shis.bosai.go.jp/wms-landslide
*
* 4301 Tokyo
* 4326 WGS84
* 4612 JGS2000
*
*Contour structure (sliding cliffs and lateral cliffs) L-V3-S100
*Contour structure (contour / boundary of moving body) (moving body arc) L-V3-S200
*Contour structure (contour / boundary of moving body) (moving body polygon) L-V3-S300
*Internal structure L-V3-S400
*Movement direction, etc. Main movement direction of the moving body L-V3-S500
*Center of gravity of mobile body L-V3-CENTER
*All the above landslide terrain L-V3-ALL
*
*/
public static BufferedImage getJShis(int zoom,int tileX,int tileY){
try{
String url=getURL(zoom,tileX,tileY);
BufferedImage img=ImageIO.read(new URL(url));
return img;
}catch(Exception e){
e.printStackTrace();
return null;
}
}
private static String getURL(int zoom,int tileX,int tileY){
String ret=baseURL+getBBOX(zoom,tileX,tileY)+"&CRS=EPSG:4612"+"&WIDTH=256&HEIGHT=256&LAYERS=L-V3-ALL&FORMAT=image/png&TRANSPARENT=TRUE";
return ret;
}
public static String getBBOX(int zoom,int tileX,int tileY){
String ret="BBOX=";
float[][] b=getTileBounds(zoom,tileX,tileY);
ret=ret+Float.toString(b[1][2])+","+Float.toString(b[0][1])+","+Float.toString(b[0][2])+","+Float.toString(b[1][1]);
return ret;
}
private static float[][] getTileBounds(int zoom,long tileX,long tileY){
float[] ll1=pixelToLonLatCoord(zoom,tileX*256,tileY*256);
float[] ll2=pixelToLonLatCoord(zoom,tileX*256+255,tileY*256+255);
return new float[][]{ll1,ll2};
}
private static float[] pixelToLonLatCoord(int zoom,long x,long y){
double lon=180.0*(x/Math.pow(2, zoom+7)-1);
double tmp0=((-Math.PI/Math.pow(2, zoom+7))*y);
double tmp1=atanh(Math.sin(L*Math.PI/180.0));
double lat=(180.0/Math.PI)*Math.asin(Math.tanh(tmp0+tmp1));
return new float[]{(float)zoom,(float)lon,(float)lat};
}
private static double atanh(double v){
return 0.5*Math.log((1.0+v)/(1.0-v));
}
}
■ CS 3D tile
■ JSHIS landslide terrain tile
■ CS landslide terrain tile
3. Other
This time, the purpose is local use, so in the code on Github, the generated tiles are stored in Sqlite. It also includes the tile browser code that uses jxmapviewer to check the generated tiles.
Reference material
In creating the CS 3D tile, I referred to the following materials.
-Terrain interpretation using CS 3D map (FOSS4G 2017 Tokyo hands-on) -[Characteristics of CS 3D drawing (Nagano Forestry Consultants Association)](http://www.rincon.or.jp/sinrinseibikasokukaringyosaiseikyogikai/CS%8E%E8%8F%87%8F%91PDF/2P9%81%6020. pdf)
Finally
@Wayama_ryousuke's entry "pix2pix refill! I tried to create a landslide topography distribution map from CS 3D map //qiita.com/wayama_ryousuke/items/8e6c2b091603a30d294f) " It seems to be very interesting, so it may be impossible in terms of resolution, but I would like to try pix2pix with the CS 3D tile and landslide terrain tile generated this time.
Recommended Posts