# 基本介绍
官网地址 (opens new window)
GeoTools是一个处理空间数据的工具包，其中包含构建一个完整的地理信息系统所需要的的全部工具类
# shpfile文件的读取和写入
读取shpfile
private static String findPgFieldName(String name, Map<String,Object> pg2ShpFieldMapping) {
String str = "";
for(Map.Entry<String, Object> entry : pg2ShpFieldMapping.entrySet()) {
if(name.equalsIgnoreCase(String.valueOf(entry.getValue()))) {
str = entry.getKey();
}
}
return str;
}
/**
* 导入 shp 到 PG
* @param path
* @return
*/
public static List<List<String>> readShpfile(String path) {
List<List<String>> dataList = new ArrayList<>();
ShapefileDataStore shpDataStore;
List<Map<String, String>> list = new ArrayList<>();
try {
shpDataStore = new ShapefileDataStore(new File(path).toURI().toURL());
//shpDataStore.setCharset(Charset.forName("UTF-8"));
shpDataStore.setCharset(Charset.forName("GB18030")); //根据用户数据的编码设定
String typeName = shpDataStore.getTypeNames()[0];
FeatureSource<SimpleFeatureType, SimpleFeature> featureSource = null;
featureSource = shpDataStore.getFeatureSource(typeName);
FeatureCollection<SimpleFeatureType, SimpleFeature> result = featureSource.getFeatures();
System.out.println(result.size());
FeatureIterator<SimpleFeature> itertor = result.features();
List<AttributeDescriptor> attributeDescriptorList = featureSource.getSchema().getAttributeDescriptors();
// 添加表头字段
dataList.add(attributeDescriptorList.stream().map(attributeDescriptor-> attributeDescriptor.getLocalName()).collect(Collectors.toList()));
while(itertor.hasNext()){
SimpleFeature feature = itertor.next();
Collection<Property> p = feature.getProperties();
Iterator<Property> it = p.iterator();
Map<String, String> map = new HashMap<>();
List<String> row = new ArrayList<>();
while(it.hasNext()) {
Property pro = it.next();
String shpFieldName = String.valueOf(pro.getName());
String shpFieldVal = String.valueOf(pro.getValue());
if("the_geom".equals(shpFieldName)) { // 要特殊处理一下
if(pro.getValue() instanceof MultiLineString) {
shpFieldVal = ((MultiLineString)(pro.getValue())).getGeometryN(0).toString();
}
else if(pro.getValue() instanceof MultiPolygon) {
shpFieldVal = ((MultiPolygon)(pro.getValue())).getGeometryN(0).toString();
}
}
row.add(shpFieldVal);
}
dataList.add(row);
}
itertor.close();
} catch (MalformedURLException e) {
e.printStackTrace();
} catch(IOException e) {
e.printStackTrace();
}
return dataList;
}
写入shpfile
public static void writeDataToShp(ShapefileDataStore ds, String type, List<Map<String, Object>> datas) {
int size = datas.size();
try {
FeatureWriter<SimpleFeatureType, SimpleFeature> writer = ds.getFeatureWriter(ds.getTypeNames()[0],
AUTO_COMMIT);
for (int i = 0; i < size; i++) {
SimpleFeature feature = writer.next();
for (Map.Entry<String, Object> item : datas.get(i).entrySet()) {
if ("the_geom".equals(item.getKey())) {
switch (type) {
case GeometryTypeConst.POINT:
feature.setAttribute("the_geom", WKTUtil.createPointByWKT(item.getValue().toString()));
break;
case GeometryTypeConst.LINESTRING:
feature.setAttribute("the_geom", WKTUtil.createLineByWKT(item.getValue().toString()));
break;
case GeometryTypeConst.POLYGON:
feature.setAttribute("the_geom", WKTUtil.createPolygonByWKT(item.getValue().toString()));
break;
default:
break;
}
} else {
feature.setAttribute(item.getKey(), item.getValue());
}
}
}
writer.write();
writer.close();
ds.dispose();
System.out.println("SHP数据写入成功");
} catch (Exception e) {
// TODO: handle exception
e.printStackTrace();
log.error("向shpfile文件写数据出错!");
}
}
public static ShapefileDataStore createShpFile(String filePath, String type,
LinkedHashMap<String, String> columns, String srid) {
File file = new File(filePath);
Map<String, Serializable> params = new HashMap<String, Serializable>();
ShapefileDataStore ds = null;
try {
params.put(ShapefileDataStoreFactory.URLP.key, file.toURI().toURL());
params.put("create spatial index", Boolean.TRUE);
ds = (ShapefileDataStore) new ShapefileDataStoreFactory().createNewDataStore(params);
// 定义图形信息和属性信息
SimpleFeatureTypeBuilder tb = new SimpleFeatureTypeBuilder();
// tb.setCRS(DefaultGeographicCRS.WGS84);
tb.setCRS(CRS.decode("EPSG:" + srid));
tb.setName("shapefile");
switch (type) {
case GeometryTypeConst.POINT:
tb.add("the_geom", Point.class);
break;
case GeometryTypeConst.LINESTRING:
tb.add("the_geom", LineString.class);
break;
case GeometryTypeConst.POLYGON:
tb.add("the_geom", Polygon.class);
break;
default:
break;
}
for (Map.Entry<String, String> item : columns.entrySet()) {
switch (item.getValue()) {
case "LONG":
tb.add(item.getKey(), Long.class);
break;
case "FLOAT":
tb.add(item.getKey(), Float.class);
break;
case "STRING":
tb.add(item.getKey(), String.class);
break;
default:
tb.add(item.getKey(), String.class);
break;
}
}
ds.createSchema(tb.buildFeatureType());
ds.setCharset(Charset.forName("UTF-8"));
} catch (IOException | FactoryException e) {
// TODO Auto-generated catch block
e.printStackTrace();
log.error("创建指定路径：" + filePath + " shp失败", e);
}
return ds;
}
# DEM坡度分析
DEM坡度算法参考 (opens new window) 计算指定像素坐标点的坡度
public float calcSlope(int cellX, int cellY, PlanarImage image) throws IOException {
DecimalFormat df = new DecimalFormat("#.0000");
final int[] dest = null;
int e = image.getTile(image.XToTileX(cellX), image.YToTileY(cellY)).getPixel(cellX, cellY, dest)[0];
int e1 = image.getTile(image.XToTileX(cellX - 1), image.YToTileY(cellY)).getPixel(cellX - 1, cellY, dest)[0];
int e2 = image.getTile(image.XToTileX(cellX), image.YToTileY(cellY - 1)).getPixel(cellX, cellY - 1, dest)[0];
int e3 = image.getTile(image.XToTileX(cellX + 1), image.YToTileY(cellY)).getPixel(cellX + 1, cellY, dest)[0];
int e4 = image.getTile(image.XToTileX(cellX), image.YToTileY(cellY + 1)).getPixel(cellX, cellY + 1, dest)[0];
int e5 = image.getTile(image.XToTileX(cellX - 1), image.YToTileY(cellY - 1)).getPixel(cellX - 1, cellY - 1,
dest)[0];
int e6 = image.getTile(image.XToTileX(cellX + 1), image.YToTileY(cellY - 1)).getPixel(cellX + 1, cellY - 1,
dest)[0];
int e7 = image.getTile(image.XToTileX(cellX + 1), image.YToTileY(cellY + 1)).getPixel(cellX + 1, cellY + 1,
dest)[0];
int e8 = image.getTile(image.XToTileX(cellX - 1), image.YToTileY(cellY + 1)).getPixel(cellX - 1, cellY + 1,
dest)[0];
double slopeWE = ((e8 + 2 * e1 + e5) - (e7 + 2 * e3 + e6)) / (8 * 2041.823085);// 东西方向坡度
double slopeNW = ((e7 + 2 * e4 + e8) - (e6 + 2 * e2 + e5)) / (8 * 2041.823085);// 南北方向坡度
double slope = 100*(Math.sqrt(Math.pow(slopeWE, 2) + Math.pow(slopeNW, 2)));
return Float.parseFloat(df.format(slope));
}
DEM数据测试用例
@Test
public void testCalcSlope() throws NoSuchAuthorityCodeException, FactoryException, IOException {
String path = "D:\\workData\\geotiff\\testTiff.tif";
String outputPath = "D:\\workData\\geotiff\\output.tif";
File file = new File(path);
// 设置tiff影像默认设置
Hints tiffHints = new Hints();
tiffHints.add(new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE));
// 默认坐标系EPSG:3857
// tiffHints.add(new Hints(Hints.DEFAULT_COORDINATE_REFERENCE_SYSTEM, CRS.decode("EPSG:4326")));
tiffHints.add(new Hints(Hints.DEFAULT_COORDINATE_REFERENCE_SYSTEM, DefaultGeographicCRS.WGS84));
GeoTiffReader reader = new GeoTiffReader(file, tiffHints);
GridCoverage2D coverage = reader.read(null);
Envelope env = coverage.getEnvelope();
PlanarImage image = (PlanarImage) coverage.getRenderedImage();
int width = image.getWidth(); // Image Width
int height = image.getHeight(); // Image Height
// 计算每个栅格的坡度
float[][] slopeData = new float[height][width];
for (int i = 1; i < height + 1; i++) {
for (int j = 1; j < width + 1; j++) {
float slope = SlopeUtil.INSTANCE.calcSlope(j, i, image);
slopeData[i - 1][j - 1] = slope;
}
}
GridCoverageFactory factory = new GridCoverageFactory();
GridCoverage2D outputCoverage = factory.create("test", slopeData, env);
GeoTiffWriter writer = new GeoTiffWriter(new File(outputPath));
writer.write(outputCoverage, null);
writer.dispose();
}
# 读取geojson
{"type":"Feature",
"crs":{"type":"name","properties":{"name":"EPSG:2380"}},“geometry”:{“type”:“MultiPolygon”,“coordinates”:[[[[646398.9535,3267941.9664],[649558.7196,3267895.3528],[649674.763,3265683.4124],[646387.8773,3265827.4858],[646398.9535,3267941.9664]]]]},“properties”:{“Id”:0}}
public static void readGeoJson(String jsonPath) throws IOException {
File file = new File(jsonPath);
FileInputStream fileInputStream = new FileInputStream(file);
// 这里可能存在问题，如果是超大文件怎么办，一次性读入会不会报内存
// 解决方案是不是对大文件进行拆分
GeoJSONReader geoJSONReader = new GeoJSONReader(fileInputStream);
SimpleFeatureCollection featureCollection = geoJSONReader.getFeatures();
SimpleFeatureIterator iterator = featureCollection.features();
List<AttributeDescriptor> attributes = featureCollection.getSchema().getAttributeDescriptors();
while (iterator.hasNext()) {
SimpleFeature simpleFeature = iterator.next();
System.out.println();
attributes.stream().forEach((a) -> {
// 依次读取这个shape中每一个属性值，当然这个属性值，可以处理其它业务
System.out.println(a.getLocalName() + ":" + simpleFeature.getAttribute(a.getLocalName()));
});
}
fileInputStream.close();
System.out.println("读取JSON完毕！");
}
# EsriJSON与GeoJSON转换工具类
package piesat.geo;
import com.alibaba.fastjson.JSON;
import com.alibaba.fastjson.JSONObject;
import org.apache.commons.lang3.StringUtils;
import java.math.BigDecimal;
import java.util.*;
public class EsriGeoJsonUtil {
public static String esri2geo(String ersiJson){
Map geoMap = new HashMap();
try {
List geoFs = new ArrayList();
geoMap.put("type", "FeatureCollection");
Map esriMap = (Map) JSON.parse(ersiJson);
Object esriFs = esriMap.get("features");
if(esriFs instanceof List){
esriFs = (List<Map<String, Object>>) esriFs;
for(int i=0; i< ((List) esriFs).size(); i++){
Map esriF = (Map) ((List) esriFs).get(i);
Map geoF = new HashMap();
geoF.put("type", "Feature");
geoF.put("properties", esriF.get("attributes"));
Map<String, Object> geometry = (Map<String, Object>) esriF.get("geometry");
if(null != geometry.get("x")){
geoF.put("geometry", geoPoint(geometry));
}else if(null != geometry.get("points")){
geoF.put("geometry", geoPoints(geometry));
}else if(null != geometry.get("paths")){
geoF.put("geometry", geoLine(geometry));
}else if(null != geometry.get("rings")){
geoF.put("geometry", geoPoly(geometry));
}
geoFs.add(geoF);
}
geoMap.put("features", geoFs);
}
}catch (Exception e){
e.printStackTrace();
}
return new JSONObject(geoMap).toString();
}
public static String geo2ersi(String geoJson, String idAttribute){
Map esriMap = new HashMap();
try {
Map geoMap = (Map) JSON.parse(geoJson);
esriMap = getEsriGeo(geoMap, idAttribute);
Map spatialReference = new HashMap();
spatialReference.put("wkid", 4326);
esriMap.put("spatialReference",spatialReference);
}catch (Exception e){
e.printStackTrace();
}
return new JSONObject(esriMap).toString();
}
public static Map getEsriGeo(Map geoMap, String idAttribute){
Map esriMap = new HashMap();
idAttribute = StringUtils.isNotEmpty(idAttribute)? idAttribute:"OBJECTID";
String type = geoMap.get("type").toString();
switch (type){
case "Point":
List<BigDecimal> coords = (List<BigDecimal>) geoMap.get("coordinates");
esriMap.put("x", coords.get(0));
esriMap.put("y", coords.get(1));
break;
case "MultiPoint":
esriMap.put("points",geoMap.get("coordinates"));
break;
case "LineString":
List<Object> coordsList = new ArrayList<>();
coordsList.add(geoMap.get("coordinates"));
esriMap.put("paths",coordsList);
break;
case "MultiLineString":
esriMap.put("paths",geoMap.get("coordinates"));
break;
case "Polygon":
List<List<List<BigDecimal>>> coordinates = (List<List<List<BigDecimal>>>) geoMap.get("coordinates");
List<List<List<BigDecimal>>> rings = orientRings(coordinates);
esriMap.put("rings", rings);
break;
case "MultiPolygon":
List<List<List<List<BigDecimal>>>> mcoordinates = (List<List<List<List<BigDecimal>>>>) geoMap.get("coordinates");
List<List<List<BigDecimal>>> mrings = flattenMultiPolygonRings(mcoordinates);
esriMap.put("rings", mrings);
break;
case "Feature":
if (null != geoMap.get("geometry")) {
Map geometry = getEsriGeo((Map) geoMap.get("geometry"), idAttribute);
esriMap.put("geometry", geometry);
}
if(null != geoMap.get("properties")){
Map properties = (Map) geoMap.get("properties");
if (null != geoMap.get("id")) {
properties.put(idAttribute, geoMap.get("id"));
}
esriMap.put("attributes", properties);
}
break;
case "FeatureCollection":
List<Object> esriFs = new ArrayList<>();
List<Map> features = (List<Map>) geoMap.get("features");
for (int i = 0; i < features.size(); i++) {
esriFs.add(getEsriGeo(features.get(i), idAttribute));
}
esriMap.put("features", esriFs);
esriMap.put("geometryType","esriGeometryPolygon");
break;
case "GeometryCollection":
List<Object> esriFsc = new ArrayList<>();
List<Map> geometries = (List<Map>) geoMap.get("geometries");
for (int i = 0; i < geometries.size(); i++) {
esriFsc.add(getEsriGeo(geometries.get(i), idAttribute));
}
esriMap.put("geometries", esriFsc);
esriMap.put("geometryType","esriGeometryPolygon");
break;
}
return esriMap;
}
public static Map geoPoint(Map<String, Object> geometry){
Map geo = new HashMap();
geo.put("type", "point");
BigDecimal x = (BigDecimal) geometry.get("x");
BigDecimal y = (BigDecimal) geometry.get("y");
List<BigDecimal> coords = new ArrayList<>();
coords.add(x);
coords.add(y);
geo.put("coordinates", coords);
return geo;
}
public static Map geoPoints(Map<String, Object> geometry){
Map geo = new HashMap();
List<Object> points = (List<Object>) geometry.get("points");
if(points.size()==1){
geo.put("type", "Point");
geo.put("coordinates", points.get(0));
}else{
geo.put("type", "MultiPoint");
geo.put("coordinates", points);
}
return geo;
}
public static Map geoLine(Map<String, Object> geometry){
Map geo = new HashMap();
List<Object> paths = (List<Object>) geometry.get("paths");
if(paths.size()==1){
geo.put("type", "LineString");
geo.put("coordinates", paths.get(0));
}else{
geo.put("type", "MultiLineString");
geo.put("coordinates", paths);
}
return geo;
}
public static Map geoPoly(Map<String, Object> geometry){
Map geo = new HashMap();
List<List<List<BigDecimal>>> rings = (List<List<List<BigDecimal>>>) geometry.get("rings");
if(rings.size()==1){
geo.put("type", "Polygon");
geo.put("coordinates", rings);
}else{
List<List<List<List<BigDecimal>>>> coords = new ArrayList();
String type = "";
int len = coords.size() - 1;
for(int i=0; i< rings.size(); i++){
if(ringIsClockwise( rings.get(i))){
List<List<List<BigDecimal>>> item = new ArrayList<>();
item.add(rings.get(i));
coords.add(item);
len++;
}else{
coords.get(len).add(rings.get(i));
}
}
if(coords.size() == 1){
type="Polygon";
}else{
type="MultiPolygon";
}
geo.put("type", type);
geo.put("coordinates",coords.size()==1?coords.get(0): coords);
}
return geo;
}
public static boolean ringIsClockwise(List<List<BigDecimal>> rings){
int total = 0;
List<BigDecimal> pt1 = null;
List<BigDecimal> pt2 = null;
for(int i=0; i< rings.size()-1; i++){
pt1 = rings.get(i);
pt2 = rings.get(i+1);
total += (pt2.get(0).doubleValue() - pt1.get(0).doubleValue())* (pt2.get(1).doubleValue() + pt1.get(1).doubleValue());
}
return total>=0;
}
public static List<List<List<BigDecimal>>> orientRings ( List<List<List<BigDecimal>>> polygon) {
List<List<List<BigDecimal>>> ringsList = new ArrayList<>();
List<List<BigDecimal>> outerRing = closeRing(polygon.get(0));
if (outerRing.size() >= 4) {
if (!ringIsClockwise(outerRing)) {
Collections.reverse(outerRing);
}
ringsList.add(outerRing);
polygon.remove(0);
for (int i = 0; i < polygon.size(); i++) {
List<List<BigDecimal>> hole = closeRing(polygon.get(i));
if (hole.size() >= 4) {
if (ringIsClockwise(hole)) {
Collections.reverse(hole);
}
ringsList.add(hole);
}
}
}
return ringsList;
}
public static List<List<BigDecimal>> closeRing (List<List<BigDecimal>> coords) {
if (!pointsEqual(coords.get(0), coords.get(coords.size()-1))) {
coords.add(coords.get(0));
}
return coords;
}
public static boolean pointsEqual (List<BigDecimal> a, List<BigDecimal> b) {
for (int i = 0; i < a.size(); i++) {
if (a.get(i).compareTo(b.get(i)) != 0) {
return false;
}
}
return true;
}
public static List<List<List<BigDecimal>>> flattenMultiPolygonRings (List<List<List<List<BigDecimal>>>> rings) {
List<List<List<BigDecimal>>> polygonList = new ArrayList<>();
for (int i = 0; i < rings.size(); i++) {
List<List<List<BigDecimal>>> polygons = orientRings(rings.get(i));
for (int x = polygons.size() - 1; x >= 0; x--) {
List<List<BigDecimal>> polygon = polygons.get(x);
polygonList.add(polygon);
}
}
return polygonList;
}
}
# geojson shp 通过geopandas 相互转换
import geopandas as gpd
# shp到geojson
data = gpd.read_file('data.shp')
data.to_file("data.json", driver='GeoJSON', encoding="utf-8")
# geojson到shp
data = gpd.read_file('data.json')
data.to_file('data.shp', driver='ESRI Shapefile', encoding='utf-8')