add UTM projection

This commit is contained in:
wyq 2024-07-20 15:32:42 +08:00
parent 2447a89fe1
commit e899d6d644
18 changed files with 221 additions and 44 deletions

View File

@ -329,9 +329,10 @@ public class MapGridLine extends GridLine {
continue;
}
if (graphic.getNumGraphics() > 1) {
points = (List<PointD>) graphic.getGraphicN(0).getShape().getPoints();
List<PointD> points1 = (List<PointD>) graphic.getGraphicN(1).getShape().getPoints();
PolylineShape polylineShape = (PolylineShape) graphic.getShape();
if (polylineShape.getPartNum() > 1) {
points = (List<PointD>) polylineShape.getPolylines().get(0).getPointList();
List<PointD> points1 = (List<PointD>) polylineShape.getPolylines().get(1).getPointList();
Collections.reverse(points1);
points.addAll(points1);
line = new PolylineZShape();
@ -417,6 +418,10 @@ public class MapGridLine extends GridLine {
if (!aGL.isLongitude()) {
aGL.setLabDirection(Direction.North);
} else {
if (aGL.getLabDirection() == Direction.North) {
continue;
}
if (aGL.getCoord().Y > 0 && Math.abs(aGL.getCoord().X) < 1000) {
continue;
}

View File

@ -940,7 +940,7 @@ public class MIMath {
* @return Values
*/
public static double[] getIntervalValues(double min, double max, int n) {
return getIntervalValues(min, max, n, false);
return getIntervalValues(min, max, n, true);
}
/**

View File

@ -1013,6 +1013,9 @@ public class GeoComputation {
}
q1 = q2;
}
if (j == borderList.size()) {
j = j - 1;
}
GridLabel aGL = new GridLabel();
aGL.setLongitude(isVertical);
aGL.setBorder(true);
@ -1063,6 +1066,9 @@ public class GeoComputation {
}
q1 = q2;
}
if (j == borderList.size()) {
j = j - 1;
}
GridLabel aGL = new GridLabel();
aGL.setBorder(true);
aGL.setLongitude(isVertical);

View File

@ -64,14 +64,15 @@ public class GeometryUtil {
return cET;
} else {
Extent cET = new Extent();
Extent cET = null;
for (int i = 0; i < PList.size(); i++) {
PointD aP = PList.get(i);
if (i == 0) {
cET.minX = aP.X;
cET.maxX = aP.X;
cET.minY = aP.Y;
cET.maxY = aP.Y;
if (Double.isInfinite(aP.X) || Double.isInfinite(aP.Y)) {
continue;
}
if (cET == null) {
cET = new Extent(aP.X, aP.X, aP.Y, aP.Y);
} else {
if (cET.minX > aP.X) {
cET.minX = aP.X;

View File

@ -1,32 +1,36 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<MeteoInfo File="milconfig.xml" Type="configurefile">
<Path OpenPath="D:\Working\MIScript\Jython\mis\chart\axis">
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\burf"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\grib"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\hdf"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\radar"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\satellite\modis"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\satellite"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\gui"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\plot"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\traj"/>
<Path OpenPath="D:\Working\MIScript\Jython\mis\map\projection">
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\chart\axes"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\chart"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\chart\axis"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\map\topology"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo\eof"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo\calc"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math\stats"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\chart"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\chart\subplot"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\scatter"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\map"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\map\projection"/>
</Path>
<File>
<OpenedFiles>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\plot\plot_cdata_3.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\traj\plot_basic.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\chart\axis\axis_color.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\map\projection\stere_proj_1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\map\projection\aeqd_proj.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\map\projection\utm_proj.py"/>
</OpenedFiles>
<RecentFiles>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\plot\plot_cdata_3.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\traj\plot_basic.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\chart\axis\axis_color.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\map\projection\stere_proj_1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\map\projection\aeqd_proj.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\map\projection\utm_proj.py"/>
</RecentFiles>
</File>
<Font>
@ -34,5 +38,5 @@
</Font>
<LookFeel DockWindowDecorated="true" LafDecorated="true" Name="FlatDarkLaf"/>
<Figure DoubleBuffering="true"/>
<Startup MainFormLocation="-7,0" MainFormSize="1396,843"/>
<Startup MainFormLocation="-7,-7" MainFormSize="1293,765"/>
</MeteoInfo>

View File

@ -14,6 +14,7 @@ from org.meteoinfo.projection.info import Wagner3 as JWagner3
from org.meteoinfo.projection.info import Airy as JAiry
from org.meteoinfo.projection.info import Aitoff as JAitoff
from org.meteoinfo.projection.info import August as JAugust
from org.meteoinfo.projection.info import UTM as JUTM
from org.locationtech.proj4j import CRSFactory
@ -21,7 +22,7 @@ __all__ = ['AlbersEqualArea','Airy','Aitoff','August','AzimuthalEquidistant','Eq
'Geostationary','Hammer','LambertAzimuthalEqualArea','LambertConformal','LambertEqualArea',
'Mercator','Mollweide','NorthPolarStereo','Orthographic',
'PlateCarree','Robinson','Sinusoidal','SouthPolarStereo','Stereographic',
'TransverseMercator','Wagner3']
'TransverseMercator','UTM','Wagner3']
crs_factory = CRSFactory()
@ -110,6 +111,33 @@ class Mercator(JMercator):
JMercator.__init__(self, crs)
class UTM(JUTM):
"""
Universal Transverse Mercator projection.
"""
def __init__(self, zone, southern_hemisphere=False):
"""
Parameters
----------
zone
The numeric zone of the UTM required.
southern_hemisphere: optional
Set to True if the zone is in the southern hemisphere. Defaults to
False.
"""
proj4_params = ['+proj=utm',
'+units=m',
'+zone=' + str(zone)]
if southern_hemisphere:
proj4_params.append('+south')
crs = crs_factory.createFromParameters('custom', proj4_params)
JUTM.__init__(self, crs)
class AlbersEqualArea(Albers):
"""
An Albers Equal Area projection.

View File

@ -212,7 +212,7 @@ def inpolygon(x, y, polygon):
r = NDArray(GeometryUtil.inPolygon(x._array, y._array, x_p._array, y_p._array))
else:
if isinstance(polygon, MILayer):
polygon = polygon.shapes()
polygon = polygon.shapes
elif isinstance(polygon, PolygonShape):
polygon = [polygon]
r = NDArray(GeometryUtil.inPolygon(x._array, y._array, polygon))

View File

@ -1333,7 +1333,7 @@ def mean(x, axis=None, keepdims=False):
for xx in x:
a.append(xx.asarray())
r = ArrayMath.mean(a)
return x0.array_wrap(r)
return x[0].array_wrap(r)
elif isinstance(x[0], PyStationData):
a = []
for xx in x:

View File

@ -852,10 +852,6 @@ def axesm(*args, **kwargs):
:param rightaxis: (*boolean*) Optional, set right axis visible or not. Default is ``True`` .
:param xyscale: (*int*) Optional, set scale of x and y axis, default is 1. It is only
valid in longlat projection.
:param gridlabel: (*boolean*) Optional, set axis tick labels visible or not. Default is ``True`` .
:param gridline: (*boolean*) Optional, set grid line visible or not. Default is ``False`` .
:param griddx: (*float*) Optional, set x grid line interval. Default is 10 degree.
:param griddy: (*float*) Optional, set y grid line interval. Default is 10 degree.
:param frameon: (*boolean*) Optional, set frame visible or not. Default is ``False`` for lon/lat
projection, ortherwise is ``True``.
:param tickfontname: (*string*) Optional, set axis tick labels font name. Default is ``Arial`` .

View File

@ -8,10 +8,7 @@ package org.meteoinfo.ndarray.math;
import org.meteoinfo.ndarray.*;
import java.time.LocalDateTime;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
import java.util.*;
/**
*
@ -46,7 +43,44 @@ public class ArrayMath {
}
}
private static DataType commonType(DataType aType, DataType bType) {
/**
* Get common data type
* @param dataTypes Data type list
* @return Common data type
*/
public static DataType commonType(List<DataType> dataTypes) {
DataType cDataType = dataTypes.get(0);
if (dataTypes.size() > 1) {
for (int i = 1; i < dataTypes.size(); i++) {
DataType dataType = dataTypes.get(i);
cDataType = commonType(cDataType, dataType);
}
}
return cDataType;
}
/**
* Get common data type
* @param arrays Array list
* @return Common data type
*/
public static DataType commonTypeArrays(List<Array> arrays) {
List<DataType> dataTypes = new ArrayList<>();
for (Array array : arrays) {
dataTypes.add(array.getDataType());
}
return commonType(dataTypes);
}
/**
* Get common data type
* @param aType Data type a
* @param bType Data type b
* @return Common data type
*/
public static DataType commonType(DataType aType, DataType bType) {
if (aType == bType) {
return aType;
}

View File

@ -2083,7 +2083,8 @@ public class ArrayUtil {
}
int[] shape = arrays.get(0).getShape();
shape[axis] = len;
Array r = Array.factory(arrays.get(0).getDataType(), shape);
DataType dataType = ArrayMath.commonTypeArrays(arrays);
Array r = Array.factory(dataType, shape);
int[] current;
IndexIterator ii = r.getIndexIterator();
Index index;
@ -2133,7 +2134,9 @@ public class ArrayUtil {
nshape[i] = shape[i];
}
}
Array r = Array.factory(a.getDataType(), nshape);
DataType dataType = ArrayMath.commonType(a.getDataType(), b.getDataType());
Array r = Array.factory(dataType, nshape);
Index indexr = r.getIndex();
Index indexa = a.getIndex();
Index indexb = b.getIndex();
@ -6017,6 +6020,4 @@ public class ArrayUtil {
}
// </editor-fold>
// <editor-fold desc="Time average">
// </editor-fold>
}

View File

@ -142,6 +142,9 @@ package org.meteoinfo.projection;
case "Transverse Mercator":
projInfo = new TransverseMercator(crs);
break;
case "Extended Transverse Mercator":
projInfo = new UTM(crs);
break;
case "Wagner III":
projInfo = new Wagner3(crs);
break;

View File

@ -38,6 +38,7 @@ public enum ProjectionNames {
Geostationary_Satellite("geos"),
Oblique_Stereographic_Alternative("sterea"),
Transverse_Mercator("tmerc"),
UTM("utm"),
Sinusoidal("sinu"),
Cylindrical_Equal_Area("cea"),
Hammer_Eckert("hammer"),

View File

@ -176,7 +176,12 @@ public class ProjectionUtil {
PointD p;
for (int i = 0; i < Y.length; i++) {
for (int j = 0; j < X.length; j++) {
p = Reproject.reprojectPoint(X[j], Y[i], fromProj, toProj);
try {
p = Reproject.reprojectPoint(X[j], Y[i], fromProj, toProj);
} catch (Exception e) {
continue;
}
if (Double.isNaN(p.X) || Double.isNaN(p.Y)) {
continue;
}
@ -744,6 +749,8 @@ public class ProjectionUtil {
}
break;
case Mercator:
case Transverse_Mercator:
case UTM:
if (pointShape.getPoint().Y > cutoff || pointShape.getPoint().Y < -cutoff) {
return null;
}
@ -792,6 +799,8 @@ public class ProjectionUtil {
}
break;
case Mercator:
case Transverse_Mercator:
case UTM:
if (lineShape.getExtent().maxY > cutoff) {
lineShape = GeoComputation.clipPolylineShape_Lat(lineShape, cutoff, false);
}
@ -872,6 +881,8 @@ public class ProjectionUtil {
}
break;
case Mercator:
case Transverse_Mercator:
case UTM:
if (aPGS.getExtent().maxY > cutoff) {
aPGS = GeoComputation.clipPolygonShape_Lat(aPGS, cutoff, false);
}

View File

@ -0,0 +1,87 @@
/* Copyright 2012 - Yaqiang Wang,
* yaqiang.wang@gmail.com
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation; either version 2.1 of the License, or (at
* your option) any later version.
*
* This library is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser
* General Public License for more details.
*/
package org.meteoinfo.projection.info;
import org.locationtech.proj4j.CoordinateReferenceSystem;
import org.meteoinfo.common.PointD;
import org.meteoinfo.geometry.shape.PolygonShape;
import org.meteoinfo.projection.ProjectionInfo;
import org.meteoinfo.projection.ProjectionNames;
import java.util.ArrayList;
import java.util.List;
/**
*
* @author Yaqiang Wang
*/
public class UTM extends ProjectionInfo {
// <editor-fold desc="Variables">
// </editor-fold>
// <editor-fold desc="Constructor">
/**
* Construction
*
* @param crs Coordinate reference system
*/
public UTM(CoordinateReferenceSystem crs) {
this.crs = crs;
this.cutoff = 85.0511f;
updateBoundary();
}
// </editor-fold>
// <editor-fold desc="Get Set Methods">
/**
* Get projection name
*
* @return Projection name
*/
@Override
public ProjectionNames getProjectionName() {
return ProjectionNames.UTM;
}
// </editor-fold>
// <editor-fold desc="Methods">
/**
* Set latitude cutoff
*
* @param value Latitude cutoff
*/
@Override
public void setCutoff(float value) {
this.cutoff = value;
this.updateBoundary();
}
@Override
public void updateBoundary() {
double x0 = -2e7;
double x1 = 2e7;
double y0 = -1e7;
double y1 = 1e7;
List<PointD> points = new ArrayList<>();
points.add(new PointD(x0, y0));
points.add(new PointD(x1, y0));
points.add(new PointD(x1, y1));
points.add(new PointD(x0, y1));
points.add(new PointD(x0, y0));
PolygonShape ps = new PolygonShape();
ps.setPoints(points);
this.boundary = ps;
}
// </editor-fold>
}