add polygonindex function

This commit is contained in:
wyq 2023-02-14 15:28:24 +08:00
parent 7fac130689
commit 16a9d12261
13 changed files with 315 additions and 130 deletions

View File

@ -67,7 +67,7 @@ import java.util.zip.ZipInputStream;
public static String getVersion(){
String version = GlobalUtil.class.getPackage().getImplementationVersion();
if (version == null || version.equals("")) {
version = "3.5.5";
version = "3.5.6";
}
return version;
}

View File

@ -1183,7 +1183,7 @@ public class VectorLayer extends MapLayer {
aField.setColumnName(aField.getColumnName() + "_1");
}
}
_attributeTable.getTable().getColumns().add(aField);
_attributeTable.getTable().addColumn(aField);
}
/**

View File

@ -243,6 +243,28 @@ public class GeoComputation {
return isIn;
}
/**
* Get polygon index from polygons by a point
* @param polygons The polygons
* @param point The point
* @return Polygon index - the point is inside the polygon
*/
public static int polygonIndex(List<PolygonShape> polygons, PointD point) {
int idx = -1;
Extent ext = GeometryUtil.getExtent(polygons);
if (MIMath.pointInExtent(point, ext)) {
PolygonShape polygonShape;
for (int i = 0; i < polygons.size(); i++) {
polygonShape = polygons.get(i);
if (pointInPolygon(polygonShape, point)) {
return i;
}
}
}
return idx;
}
// /**
// * Determine if a point loacted in a polygon layer
// *

View File

@ -447,20 +447,27 @@ public class GeometryUtil {
return inPolygon(x, y, shapes);
}
// /**
// * Maskout function
// *
// * @param a Array a
// * @param x X dimension values
// * @param y Y dimension values
// * @param layer VectorLayer
// * @param missingValue Missing value
// * @return Result array with cell values of missing outside polygons
// */
// public static Array maskout(Array a, List<Number> x, List<Number> y, VectorLayer layer, Number missingValue) {
// List<PolygonShape> polygons = (List<PolygonShape>) layer.getShapes();
// return maskout(a, x, y, polygons, missingValue);
// }
/**
* Find polygon index for each point
* @param x X coordinate of the points
* @param y Y coordinate of the points
* @param polygons The polygons
* @return Polygon index for each point
*/
public static Array polygonIndex(Array x, Array y, List<PolygonShape> polygons) {
Array r = Array.factory(DataType.INT, x.getShape());
IndexIterator xIter = x.getIndexIterator();
IndexIterator yIter = y.getIndexIterator();
IndexIterator rIter = r.getIndexIterator();
int idx;
while (rIter.hasNext()){
idx = GeoComputation.polygonIndex(polygons, new PointD(xIter.getDoubleNext(),
yIter.getDoubleNext()));
rIter.setIntNext(idx);
}
return r;
}
/**
* Maskout function

View File

@ -1,14 +1,7 @@
<?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\bar">
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\radar"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\scatter"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\quiver"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\surf"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math\interpolate"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\isosurface"/>
<Path OpenPath="D:\Working\MIScript\Jython\mis\dataframe">
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\particles"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\wind"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\3d_earth"/>
@ -16,19 +9,24 @@
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\contour"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\bar"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\map\geod"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\map\geoshow"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\map\topology"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\map"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\map\maskout"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\dataframe"/>
</Path>
<File>
<OpenedFiles>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\surf\surf_pumpkin_1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\3d_earth\scatter_sphere.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\integrate\odeint_1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\integrate\odeint_lorenz.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\dataframe\dataframe_groupby_1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\map\maskout\polygonindex_1.py"/>
</OpenedFiles>
<RecentFiles>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\surf\surf_pumpkin_1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\3d_earth\scatter_sphere.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\common_math\integrate\odeint_1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\common_math\integrate\odeint_lorenz.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\dataframe\dataframe_groupby_1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\map\maskout\polygonindex_1.py"/>
</RecentFiles>
</File>
<Font>
@ -36,5 +34,5 @@
</Font>
<LookFeel DockWindowDecorated="true" LafDecorated="true" Name="FlatDarkLaf"/>
<Figure DoubleBuffering="true"/>
<Startup MainFormLocation="-7,0" MainFormSize="1367,792"/>
<Startup MainFormLocation="-7,0" MainFormSize="1381,848"/>
</MeteoInfo>

View File

@ -33,9 +33,10 @@ from java.util import ArrayList
__all__ = [
'arrayinpolygon', 'bilwrite', 'circle', 'convert_encoding_dbf', 'distance', 'georead', 'geotiffread', 'gridarea',
'maplayer', 'inpolygon', 'maskin', 'maskout', 'polyarea', 'polygon', 'rmaskin', 'rmaskout', 'shaperead',
'projinfo','project','projectxy','reproject','reproject_image'
'polygonindex', 'projinfo', 'project', 'projectxy', 'reproject', 'reproject_image'
]
def shaperead(fn, encoding=None):
"""
Returns a layer readed from a shape file.
@ -65,9 +66,11 @@ def shaperead(fn, encoding=None):
except:
raise
else:
print 'File not exists: ' + fn
print
'File not exists: ' + fn
raise
def georead(fn, encoding=None):
"""
Returns a layer read from a supported geo-data file.
@ -100,6 +103,7 @@ def georead(fn, encoding=None):
print('File not exists: ' + fn)
raise IOError
def geotiffread(filename):
"""
Return data array from a GeoTiff data file.
@ -113,6 +117,7 @@ def geotiffread(filename):
r = geotiff.readArray()
return NDArray(r)
def convert_encoding_dbf(filename, fromencoding, toencoding):
"""
Convert encoding of a dBase file (.dbf).
@ -131,6 +136,7 @@ def convert_encoding_dbf(filename, fromencoding, toencoding):
atable.setEncoding(toencoding)
atable.save()
def maplayer(shapetype='polygon'):
"""
Create a new map layer.
@ -142,6 +148,7 @@ def maplayer(shapetype='polygon'):
"""
return MILayer(shapetype=shapetype)
def polygon(x, y=None):
"""
Create polygon from coordinate data.
@ -161,6 +168,7 @@ def polygon(x, y=None):
polygon = ShapeUtil.createPolygonShape(x, y)
return polygon
def circle(xy, radius=5):
"""
Create a circle patch
@ -173,6 +181,7 @@ def circle(xy, radius=5):
cc = ShapeUtil.createCircleShape(xy[0], xy[1], radius)
return cc
def inpolygon(x, y, polygon):
"""
Check if x/y points are inside a polygon or not.
@ -211,6 +220,7 @@ def inpolygon(x, y, polygon):
else:
return r
def arrayinpolygon(a, polygon, x=None, y=None):
"""
Set array element value as 1 if inside a polygon or set value as -1.
@ -226,6 +236,7 @@ def arrayinpolygon(a, polygon, x=None, y=None):
if x is None or y is None:
x = self.dimvalue(1)
y = self.dimvalue(0)
if not x is None and not y is None:
if isinstance(polygon, tuple):
x_p = polygon[0]
@ -242,6 +253,27 @@ def arrayinpolygon(a, polygon, x=None, y=None):
else:
return None
def polygonindex(x, y, polygons):
"""
Find polygon index by x, y coordinates.
:param x: (*float*) X coordinate of the point.
:param y: (*float*) Y coordinate of the point.
:param polygons: (*list of PolygonShape*) The polygons
:returns: (*array_like*) Result array.
"""
if isinstance(x, (tuple, list)):
x = np.array(x)
if isinstance(x, (tuple, list)):
y = np.array(y)
r = GeometryUtil.polygonIndex(x.jarray, y.jarray, polygons)
return NDArray(r)
def distance(*args, **kwargs):
"""
Get distance of a line.
@ -264,6 +296,7 @@ def distance(*args, **kwargs):
r = GeoComputation.getDistance(x, y, islonlat)
return r
def polyarea(*args, **kwargs):
"""
Calculate area of polygon.
@ -288,6 +321,7 @@ def polyarea(*args, **kwargs):
r = GeoComputation.getArea(x, y, islonlat)
return r
def gridarea(x_orig, x_cell, x_num, y_orig, y_cell, y_num, islonlat=False,
allcell=True, earth_radius=None):
"""
@ -312,6 +346,7 @@ def gridarea(x_orig, x_cell, x_num, y_orig, y_cell, y_num, islonlat=False,
islonlat, allcell, earth_radius)
return NDArray(a)
def maskout(data, mask, x=None, y=None):
"""
Maskout data by polygons - NaN values of elements outside polygons.
@ -357,6 +392,7 @@ def maskout(data, mask, x=None, y=None):
else:
return NDArray(r)
def rmaskout(data, x, y, mask):
"""
Maskout data by polygons - the elements outside polygons will be removed
@ -373,6 +409,7 @@ def rmaskout(data, x, y, mask):
r = GeometryUtil.maskout_Remove(data.asarray(), x.asarray(), y.asarray(), mask)
return NDArray(r[0]), NDArray(r[1]), NDArray(r[2])
def maskin(data, mask, x=None, y=None):
"""
Maskin data by polygons - NaN values of elements inside polygons.
@ -411,6 +448,7 @@ def maskin(data, mask, x=None, y=None):
else:
return NDArray(r)
def rmaskin(data, x, y, mask):
"""
Maskin data by polygons - the elements inside polygons will be removed
@ -427,6 +465,7 @@ def rmaskin(data, x, y, mask):
r = GeometryUtil.maskin_Remove(data._array, x._array, y._array, mask)
return NDArray(r[0]), NDArray(r[1]), NDArray(r[2])
def projinfo(proj4string=None, proj='longlat', **kwargs):
"""
Create a projection object with Proj.4 parameters (http://proj4.org/)
@ -525,7 +564,9 @@ def projinfo(proj4string=None, proj='longlat', **kwargs):
return ProjectionInfo.factory(projstr)
def project(x, y, fromproj=KnownCoordinateSystems.geographic.world.WGS1984, toproj=KnownCoordinateSystems.geographic.world.WGS1984):
def project(x, y, fromproj=KnownCoordinateSystems.geographic.world.WGS1984,
toproj=KnownCoordinateSystems.geographic.world.WGS1984):
"""
Project geographic coordinates from one projection to another.
@ -552,6 +593,7 @@ def project(x, y, fromproj=KnownCoordinateSystems.geographic.world.WGS1984, topr
outpt = Reproject.reprojectPoint(inpt, fromproj, toproj)
return outpt.X, outpt.Y
def projectxy(lon, lat, xnum, ynum, dx, dy, toproj, fromproj=None, pos='lowerleft'):
"""
Get projected x, y coordinates by projection and a given lon, lat coordinate.
@ -581,6 +623,7 @@ def projectxy(lon, lat, xnum, ynum, dx, dy, toproj, fromproj=None, pos='lowerlef
yy = np.arange1(lly, ynum, dy)
return xx, yy
def reproject(a, x=None, y=None, fromproj=None, xp=None, yp=None, toproj=None, method='bilinear'):
"""
Project array
@ -629,6 +672,7 @@ def reproject(a, x=None, y=None, fromproj=None, xp=None, yp=None, toproj=None, m
r = Reproject.reproject(a.asarray(), x.aslist(), y.aslist(), xp, yp, fromproj, toproj, method)
return NDArray(r)
def reproject_image(a, x=None, y=None, fromproj=None, xp=None, yp=None, toproj=None):
"""
Project image data array
@ -668,6 +712,7 @@ def reproject_image(a, x=None, y=None, fromproj=None, xp=None, yp=None, toproj=N
r = Reproject.reprojectImage(a.asarray(), x.asarray(), y.asarray(), xp, yp, fromproj, toproj)
return NDArray(r)
def bilwrite(fn, data, x, y, proj=projinfo()):
"""
Write a bil file from a 2D array.

View File

@ -18,6 +18,7 @@ from org.meteoinfo.geometry.shape import PolygonShape, ShapeTypes
from org.meteoinfo.geo.analysis import GeometryUtil
from org.meteoinfo.geo.util import GeoProjectionUtil
class MILayer(object):
"""
Map layer
@ -26,6 +27,7 @@ class MILayer(object):
:param shapetype: (*ShapeTypes*) Shape type ['point' | 'point_z' | 'line' | 'line_z' | 'polygon'
| 'polygon_z']
"""
def __init__(self, layer=None, shapetype=None):
if layer is None:
if shapetype is None:
@ -121,6 +123,16 @@ class MILayer(object):
"""
return self._layer.getAttributeTable().getEncoding()
@property
def datatable(self):
"""
Get attribute table.
:return: Attribute table.
"""
r = self._layer.getAttributeTable().getTable()
return np.datatable(r)
def gettable(self):
"""
Get attribute table.
@ -152,10 +164,25 @@ class MILayer(object):
:param fieldname: (*string*) Field name.
:param shapeindex: (*int*) Shape index.
:param value: (*object*) Cell value to be asigned.
:param value: (*object*) Cell value to be assigned.
"""
self._layer.editCellValue(fieldname, shapeindex, value)
def setfieldvalue(self, fieldname, value, index=None):
"""
Set field value.
:param fieldname: (*str*) The field name.
:param value: (*array*) The field data array.
:param index: (*array*) Optional. Field data index. Default is `None`.
"""
value = np.asarray(value)
if index is None:
self._layer.getAttributeTable().getTable().setColumnData(fieldname, value.jarray)
else:
index = np.asarray(index)
self._layer.getAttributeTable().getTable().setColumnData(fieldname, value.jarray, index.jarray)
def shapes(self):
"""
Get shapes.
@ -575,4 +602,3 @@ class MIXYListData():
self.data.addSeries(key, xdata, ydata)
else:
self.data.addSeries(key, xdata.asarray(), ydata.asarray())

View File

@ -37,6 +37,10 @@ class NDArray(object):
for i in range(1, self.ndim):
self.sizestr = self.sizestr + '*%s' % self.shape[i]
@property
def jarray(self):
return self._array
#---- shape property
def get_shape(self):
return self._shape
@ -1149,12 +1153,12 @@ class NDArray(object):
def ravel(self):
"""
Return a copy of the array collapsed into one dimension.
Return a contiguous flattened array.
:returns: (*NDArray*) A copy of the input array, flattened to one dimension.
:returns: (*NDArray*) A contiguous flattened array.
"""
shape = [self.size]
r = NDArray(self._array.reshape(shape))
r = NDArray(self._array.reshapeNoCopy(shape))
return r
def repeat(self, repeats, axis=None):

View File

@ -32,6 +32,8 @@
*/
package org.meteoinfo.ndarray;
import org.meteoinfo.common.util.JDateUtil;
import java.nio.ByteBuffer;
import java.nio.DoubleBuffer;
import java.time.LocalDateTime;
@ -181,11 +183,11 @@ public class ArrayDate extends Array {
}
public double getDouble(Index i) {
throw new ForbiddenConversionException();
return JDateUtil.toOADate(storage[i.currentElement()]);
}
public void setDouble(Index i, double value) {
throw new ForbiddenConversionException();
storage[i.currentElement()] = JDateUtil.fromOADate(value);
}
public float getFloat(Index i) {
@ -288,11 +290,11 @@ public class ArrayDate extends Array {
// trusted, assumes that individual dimension lengths have been checked
public double getDouble(int index) {
throw new ForbiddenConversionException();
return JDateUtil.toOADate(storage[index]);
}
public void setDouble(int index, double value) {
throw new ForbiddenConversionException();
storage[index] = JDateUtil.fromOADate(value);
}
public float getFloat(int index) {

View File

@ -121,7 +121,7 @@ public class DataRow {
}
/**
* Set a vlaue
* Set a value
*
* @param column The data column
* @param value The value

View File

@ -14,6 +14,7 @@
package org.meteoinfo.table;
import org.meteoinfo.common.MIMath;
import org.meteoinfo.ndarray.Array;
import org.meteoinfo.ndarray.DataType;
import org.meteoinfo.ndarray.Range;
import org.meteoinfo.table.util.TableUtil;
@ -608,6 +609,86 @@ public class DataTable {
setColumnData(col, colData);
}
/**
* Add column data
*
* @param col The column
* @param colData The column data
* @throws Exception
*/
public void setColumnData(DataColumn col, Array colData) throws Exception {
colData = colData.copyIfView();
if (this.getRowCount() == 0) {
for (int i = 0; i < colData.getSize(); i++){
DataRow row = this.addRow();
row.setValue(col, colData.getObject(i));
}
} else {
int i = 0;
for (DataRow row : this.rows) {
if (i < colData.getSize()) {
row.setValue(col, colData.getObject(i));
}
i++;
}
}
}
/**
* Add column data
*
* @param colName Column name
* @param colData The column data
* @throws Exception
*/
public void setColumnData(String colName, Array colData) throws Exception {
DataColumn col = this.findColumn(colName);
if (col == null){
System.out.println("The column not exists: " + colName + "!");
return;
}
setColumnData(col, colData);
}
/**
* Add column data
*
* @param col The column
* @param colData The column data
* @param index Data index array
* @throws Exception
*/
public void setColumnData(DataColumn col, Array colData, Array index) throws Exception {
colData = colData.copyIfView();
index = index.copyIfView();
int idx;
for (int i = 0; i < index.getSize(); i++) {
idx = index.getInt(i);
if (idx >=0 && idx < this.getRowCount())
this.rows.get(idx).setValue(col, colData.getObject(i));
}
}
/**
* Add column data
*
* @param colName Column name
* @param colData The column data
* @param index Data index array
* @throws Exception
*/
public void setColumnData(String colName, Array colData, Array index) throws Exception {
DataColumn col = this.findColumn(colName);
if (col == null){
System.out.println("The column not exists: " + colName + "!");
return;
}
setColumnData(col, colData, index);
}
/**
* Add column data
*