update graphic projection functions

This commit is contained in:
wyq 2022-02-06 16:02:21 +08:00
parent 5426cf353e
commit 913c9696bb
4 changed files with 461 additions and 27 deletions

View File

@ -1213,4 +1213,54 @@ public class MIMath {
r[1] += y;
return r;
}
/**
* Get wind direction and wind speed from U/V
*
* @param u U component
* @param v V component
* @return Wind direction and wind speed
*/
public static double[] uv2ds(double u, double v) {
double ws = Math.sqrt(u * u + v * v);
double wd;
if (ws == 0) {
wd = 0;
} else {
wd = Math.asin(u / ws) * 180 / Math.PI;
if (u <= 0 && v < 0) {
wd = 180.0 - wd;
} else if (u > 0 && v < 0) {
wd = 180.0 - wd;
} else if (u < 0 && v > 0) {
wd = 360.0 + wd;
}
wd += 180;
if (wd >= 360) {
wd -= 360;
}
}
return new double[]{wd, ws};
}
/**
* Get wind U/V components from wind direction and speed
*
* @param windDir Wind direction
* @param windSpeed Wind speed
* @return Wind U/V components
*/
public static double[] ds2uv(double windDir, double windSpeed) {
double dir;
dir = windDir + 180;
if (dir > 360) {
dir = dir - 360;
}
dir = dir * Math.PI / 180;
double u = windSpeed * Math.sin(dir);
double v = windSpeed * Math.cos(dir);
return new double[]{u, v};
}
}

View File

@ -13,9 +13,7 @@
*/
package org.meteoinfo.geometry.geoprocess;
import org.meteoinfo.common.Extent;
import org.meteoinfo.common.MIMath;
import org.meteoinfo.common.PointD;
import org.meteoinfo.common.*;
import org.meteoinfo.geometry.shape.*;
//import org.meteoinfo.map.GridLabel;
//import org.meteoinfo.math.meteo.MeteoMath;
@ -860,6 +858,245 @@ public class GeoComputation {
// return newLayer;
// }
/**
* Get grid labels of a polyline
*
* @param inPolyLine Polyline
* @param clipExtent Clipping object
* @param isVertical If is vertical
* @return Clip points
*/
public static List<GridLabel> getGridLabels(Polyline inPolyLine, Extent clipExtent, boolean isVertical) {
List<GridLabel> gridLabels = new ArrayList<>();
List<PointD> aPList = (List<PointD>) inPolyLine.getPointList();
if (!isExtentCross(inPolyLine.getExtent(), clipExtent)) {
return gridLabels;
}
int i, j;
//Judge if all points of the polyline are in the cut polygon - outline
List<List<PointD>> newLines = new ArrayList<>();
PointD p1, p2;
boolean isReversed = false;
if (pointInClipObj(clipExtent, aPList.get(0))) {
boolean isAllIn = true;
int notInIdx = 0;
for (i = 0; i < aPList.size(); i++) {
if (!pointInClipObj(clipExtent, aPList.get(i))) {
notInIdx = i;
isAllIn = false;
break;
}
}
if (!isAllIn) //Put start point outside of the cut polygon
{
if (inPolyLine.isClosed()) {
List<PointD> bPList = new ArrayList<>();
bPList.addAll(aPList.subList(notInIdx, aPList.size() - 1));
bPList.addAll(aPList.subList(1, notInIdx));
bPList.add(bPList.get(0));
newLines.add(bPList);
} else {
Collections.reverse(aPList);
newLines.add(aPList);
isReversed = true;
}
} else { //the input polygon is inside the cut polygon
p1 = aPList.get(0);
if (aPList.size() == 2)
p2 = aPList.get(1);
else
p2 = aPList.get(2);
GridLabel aGL = new GridLabel();
aGL.setLongitude(isVertical);
aGL.setBorder(false);
aGL.setCoord(p1);
if (isVertical) {
aGL.setLabDirection(Direction.South);
} else {
aGL.setLabDirection(Direction.Weast);
}
aGL.setAnge((float) MIMath.uv2ds(p2.X - p1.X, p2.Y - p1.Y)[0]);
gridLabels.add(aGL);
p1 = aPList.get(aPList.size() - 1);
if (aPList.size() == 2)
p2 = aPList.get(aPList.size() - 2);
else
p2 = aPList.get(aPList.size() - 3);
aGL = new GridLabel();
aGL.setLongitude(isVertical);
aGL.setBorder(false);
aGL.setCoord(p1);
if (isVertical) {
aGL.setLabDirection(Direction.North);
} else {
aGL.setLabDirection(Direction.East);
}
aGL.setAnge((float) MIMath.uv2ds(p2.X - p1.X, p2.Y - p1.Y)[0]);
gridLabels.add(aGL);
return gridLabels;
}
} else {
newLines.add(aPList);
}
//Prepare border point list
List<BorderPoint> borderList = new ArrayList<>();
BorderPoint aBP;
List<PointD> clipPList = getClipPointList(clipExtent);
for (PointD aP : clipPList) {
aBP = new BorderPoint();
aBP.Point = aP;
aBP.Id = -1;
borderList.add(aBP);
}
//Cutting
for (int l = 0; l < newLines.size(); l++) {
aPList = newLines.get(l);
boolean isInPolygon = pointInClipObj(clipExtent, aPList.get(0));
PointD q1, q2, IPoint = new PointD();
Line lineA, lineB;
List<PointD> newPlist = new ArrayList<>();
//Polyline bLine = new Polyline();
p1 = aPList.get(0);
int inIdx = -1, outIdx = -1;
//bool newLine = true;
int a1 = 0;
for (i = 1; i < aPList.size(); i++) {
p2 = aPList.get(i);
if (pointInClipObj(clipExtent, p2)) {
if (!isInPolygon) {
lineA = new Line();
lineA.P1 = p1;
lineA.P2 = p2;
q1 = borderList.get(0).Point;
for (j = 1; j < borderList.size(); j++) {
q2 = borderList.get(j).Point;
lineB = new Line();
lineB.P1 = q1;
lineB.P2 = q2;
if (isLineSegmentCross(lineA, lineB)) {
IPoint = getCrossPoint(lineA, lineB);
inIdx = j;
break;
}
q1 = q2;
}
GridLabel aGL = new GridLabel();
aGL.setLongitude(isVertical);
aGL.setBorder(true);
aGL.setCoord(IPoint);
if (MIMath.doubleEquals(q1.X, borderList.get(j).Point.X)) {
if (MIMath.doubleEquals(q1.X, clipExtent.minX)) {
aGL.setLabDirection(Direction.Weast);
} else {
aGL.setLabDirection(Direction.East);
}
} else {
if (MIMath.doubleEquals(q1.Y, clipExtent.minY)) {
aGL.setLabDirection(Direction.South);
} else {
aGL.setLabDirection(Direction.North);
}
}
if (isVertical) {
if (aGL.getLabDirection() == Direction.South || aGL.getLabDirection() == Direction.North) {
gridLabels.add(aGL);
}
} else {
if (aGL.getLabDirection() == Direction.East || aGL.getLabDirection() == Direction.Weast) {
gridLabels.add(aGL);
}
}
}
newPlist.add(aPList.get(i));
isInPolygon = true;
} else {
if (isInPolygon) {
lineA = new Line();
lineA.P1 = p1;
lineA.P2 = p2;
q1 = borderList.get(0).Point;
for (j = 1; j < borderList.size(); j++) {
q2 = borderList.get(j).Point;
lineB = new Line();
lineB.P1 = q1;
lineB.P2 = q2;
if (isLineSegmentCross(lineA, lineB)) {
IPoint = getCrossPoint(lineA, lineB);
outIdx = j;
a1 = inIdx;
break;
}
q1 = q2;
}
GridLabel aGL = new GridLabel();
aGL.setBorder(true);
aGL.setLongitude(isVertical);
aGL.setCoord(IPoint);
if (MIMath.doubleEquals(q1.X, borderList.get(j).Point.X)) {
if (MIMath.doubleEquals(q1.X, clipExtent.minX)) {
aGL.setLabDirection(Direction.Weast);
} else {
aGL.setLabDirection(Direction.East);
}
} else {
if (MIMath.doubleEquals(q1.Y, clipExtent.minY)) {
aGL.setLabDirection(Direction.South);
} else {
aGL.setLabDirection(Direction.North);
}
}
if (isVertical) {
if (aGL.getLabDirection() == Direction.South || aGL.getLabDirection() == Direction.North) {
gridLabels.add(aGL);
}
} else {
if (aGL.getLabDirection() == Direction.East || aGL.getLabDirection() == Direction.Weast) {
gridLabels.add(aGL);
}
}
isInPolygon = false;
newPlist = new ArrayList<>();
}
}
p1 = p2;
}
if (isInPolygon && newPlist.size() > 1) {
GridLabel aGL = new GridLabel();
aGL.setLongitude(isVertical);
aGL.setBorder(false);
aGL.setCoord(newPlist.get(newPlist.size() - 1));
if (isVertical) {
if (isReversed) {
aGL.setLabDirection(Direction.South);
} else {
aGL.setLabDirection(Direction.North);
}
} else {
if (isReversed) {
aGL.setLabDirection(Direction.Weast);
} else {
aGL.setLabDirection(Direction.East);
}
}
gridLabels.add(aGL);
}
}
return gridLabels;
}
/**
* Clip a shape
*

View File

@ -1,10 +1,6 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<MeteoInfo File="milconfig.xml" Type="configurefile">
<Path OpenPath="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\volume">
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo\wrf"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\chart"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\scatter"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\arrow"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\wind"/>
@ -15,18 +11,28 @@
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d_earth"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\geoshow"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\chart"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\map"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\map\projection"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\volume"/>
</Path>
<File>
<OpenedFiles>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\volume\volumeplot_2.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\volume\volumeplot_4.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\geoshow\typhoon_map_volume.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\geoshow\geoshow_proj.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\volume\plot_cuace_3d_volume_opacity.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\volume\volumeplot_sagittal_specular.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\map\projection\lcc_proj_image.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\volume\volumeplot_1.py"/>
</OpenedFiles>
<RecentFiles>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\volume\volumeplot_2.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\volume\volumeplot_4.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\geoshow\typhoon_map_volume.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\geoshow\geoshow_proj.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\volume\plot_cuace_3d_volume_opacity.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\volume\volumeplot_sagittal_specular.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\map\projection\lcc_proj_image.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\volume\volumeplot_1.py"/>
</RecentFiles>
</File>
<Font>
@ -34,5 +40,5 @@
</Font>
<LookFeel DockWindowDecorated="true" LafDecorated="true" Name="FlatDarkLaf"/>
<Figure DoubleBuffering="true"/>
<Startup MainFormLocation="-7,-7" MainFormSize="1293,685"/>
<Startup MainFormLocation="-7,0" MainFormSize="1344,758"/>
</MeteoInfo>

View File

@ -26,6 +26,8 @@ import org.meteoinfo.ndarray.math.ArrayUtil;
import java.nio.DoubleBuffer;
import java.util.ArrayList;
import java.util.List;
import java.util.logging.Level;
import java.util.logging.Logger;
/**
*
@ -547,10 +549,12 @@ public class ProjectionUtil {
Reproject.reprojectPoints(points, fromProj, toProj, 0, points.length);
if (!Double.isNaN(points[0][0]) || !Double.isInfinite(points[0][0]) ||
!Double.isNaN(points[0][1]) || !Double.isInfinite(points[0][1])) {
//wPoint = new PointZ();
wPoint.X = points[0][0];
wPoint.Y = points[0][1];
newPoints.add(wPoint);
PointZ nPoint = new PointZ();
nPoint.X = points[0][0];
nPoint.Y = points[0][1];
nPoint.M = wPoint.M;
nPoint.Z = wPoint.Z;
newPoints.add(nPoint);
}
} catch (Exception e) {
break;
@ -638,7 +642,132 @@ public class ProjectionUtil {
}
/**
* Project polygon shape - clip the polygon is necessary
* Project point shape - clip the point shape when necessary
*
* @param pointShape A point shape
* @param fromProj From projection
* @param toProj To projection
* @return Projected point shape
*/
public static PointShape projectClipPointShape(PointShape pointShape, ProjectionInfo fromProj, ProjectionInfo toProj) {
double refLon = toProj.getCoordinateReferenceSystem().getProjection().getProjectionLongitudeDegrees();
refLon += 180;
if (refLon > 180) {
refLon = refLon - 360;
} else if (refLon < -180) {
refLon = refLon + 360;
}
float cutoff = toProj.getCutoff();
if (fromProj.getProjectionName() == ProjectionNames.LongLat) {
switch (toProj.getProjectionName()) {
case Lambert_Conformal_Conic:
if (pointShape.getPoint().Y < cutoff) {
return null;
}
break;
case North_Polar_Stereographic_Azimuthal:
if (pointShape.getPoint().Y < cutoff) {
return null;
}
break;
case South_Polar_Stereographic_Azimuthal:
if (pointShape.getPoint().Y > cutoff) {
return null;
}
break;
case Mercator:
if (pointShape.getPoint().Y > cutoff || pointShape.getPoint().Y < -cutoff) {
return null;
}
break;
}
}
pointShape = ProjectionUtil.projectPointShape(pointShape, fromProj, toProj);
return pointShape;
}
/**
* Project polyline shape - clip the polyline shape when necessary
*
* @param lineShape A polyline shape
* @param fromProj From projection
* @param toProj To projection
* @return Projected polyline shape
*/
public static List<PolylineShape> projectClipPolylineShape(PolylineShape lineShape, ProjectionInfo fromProj, ProjectionInfo toProj) {
double refLon = toProj.getCoordinateReferenceSystem().getProjection().getProjectionLongitudeDegrees();
refLon += 180;
if (refLon > 180) {
refLon = refLon - 360;
} else if (refLon < -180) {
refLon = refLon + 360;
}
float cutoff = toProj.getCutoff();
List<PolylineShape> lineShapes = new ArrayList<>();
if (fromProj.getProjectionName() == ProjectionNames.LongLat) {
switch (toProj.getProjectionName()) {
case Lambert_Conformal_Conic:
if (lineShape.getExtent().minY < cutoff) {
lineShape = GeoComputation.clipPolylineShape_Lat(lineShape, cutoff, true);
}
break;
case North_Polar_Stereographic_Azimuthal:
if (lineShape.getExtent().minY < cutoff) {
lineShape = GeoComputation.clipPolylineShape_Lat(lineShape, cutoff, true);
}
break;
case South_Polar_Stereographic_Azimuthal:
if (lineShape.getExtent().maxY > cutoff) {
lineShape = GeoComputation.clipPolylineShape_Lat(lineShape, cutoff, false);
}
break;
case Mercator:
if (lineShape.getExtent().maxY > cutoff) {
lineShape = GeoComputation.clipPolylineShape_Lat(lineShape, cutoff, false);
}
if (lineShape.getExtent().minY < -cutoff) {
lineShape = GeoComputation.clipPolylineShape_Lat(lineShape, -cutoff, true);
}
break;
}
if (lineShape == null) {
return null;
}
if (lineShape.getExtent().minX <= refLon && lineShape.getExtent().maxX >= refLon) {
switch (toProj.getProjectionName()) {
case North_Polar_Stereographic_Azimuthal:
case South_Polar_Stereographic_Azimuthal:
lineShapes.add(lineShape);
break;
default:
lineShapes.add(GeoComputation.clipPolylineShape_Lon(lineShape, refLon));
break;
}
} else {
lineShapes.add(lineShape);
}
} else {
lineShapes.add(lineShape);
}
List<PolylineShape> newPolylines = new ArrayList<>();
for (int i = 0; i < lineShapes.size(); i++) {
lineShape = lineShapes.get(i);
lineShape = ProjectionUtil.projectPolylineShape(lineShape, fromProj, toProj);
if (lineShape != null) {
newPolylines.add(lineShape);
}
}
return newPolylines;
}
/**
* Project polygon shape - clip the polygon shape when necessary
*
* @param aPGS A polygon shape
* @param fromProj From projection
@ -752,16 +881,21 @@ public class ProjectionUtil {
*/
public static Graphic projectClipGraphic(Graphic graphic, ProjectionInfo fromProj, ProjectionInfo toProj) {
if (graphic instanceof GraphicCollection) {
GraphicCollection newGCollection = new GraphicCollection();
for (Graphic aGraphic : ((GraphicCollection) graphic).getGraphics()) {
List<? extends Shape> shapes = projectClipShape(aGraphic.getShape(), fromProj, toProj);
if (shapes != null) {
aGraphic.setShape(shapes.get(0));
newGCollection.add(aGraphic);
try {
Graphic newGCollection = graphic.getClass().getDeclaredConstructor().newInstance();
for (Graphic aGraphic : ((GraphicCollection) graphic).getGraphics()) {
List<? extends Shape> shapes = projectClipShape(aGraphic.getShape(), fromProj, toProj);
if (shapes != null && shapes.size() > 0) {
aGraphic.setShape(shapes.get(0));
((GraphicCollection) newGCollection).add(aGraphic);
}
}
}
return newGCollection;
return newGCollection;
} catch (Exception ex) {
ex.printStackTrace();
return null;
}
} else {
List<? extends Shape> shapes = projectClipShape(graphic.getShape(), fromProj, toProj);
if (shapes != null) {
@ -838,11 +972,18 @@ public class ProjectionUtil {
switch (shape.getShapeType()) {
case POINT:
case POINT_M:
//shapes = projectPointShape((PointShape) shape, fromProj, toProj);
case POINT_Z:
PointShape pointShape = projectClipPointShape((PointShape) shape, fromProj, toProj);
if (pointShape != null) {
List<PointShape> pointShapes = new ArrayList<PointShape>();
pointShapes.add(pointShape);
shapes = pointShapes;
}
break;
case POLYLINE:
case POLYLINE_M:
//shapes = projectPolylineShape((PolylineShape) shape, fromProj, toProj);
case POLYLINE_Z:
shapes = projectClipPolylineShape((PolylineShape) shape, fromProj, toProj);
break;
case CURVE_LINE:
//shapes = projectCurvelineShape((CurveLineShape) shape, fromProj, toProj);