Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added sphere to plate function for meshes #944

Draft
wants to merge 1 commit into
base: master
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions __init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -234,6 +234,7 @@ def draw(self, context):
self.layout.operator("tesselation.voronoi", icon_value=icons_dict["voronoi"].icon_id, text='Voronoi')
if EARTH_SPHERE:
self.layout.operator("earth.sphere", icon="WORLD", text='lonlat to sphere')
self.layout.operator("earth.plane", icon="VIEW_PERSPECTIVE", text='sphere to lonlat')
#self.layout.operator("earth.curvature", icon="SPHERECURVE", text='Earth curvature correction')
self.layout.operator("earth.curvature", icon_value=icons_dict["curve"].icon_id, text='Earth curvature correction')

Expand Down
10 changes: 9 additions & 1 deletion core/georaster/npimg.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@

from .georef import GeoRef
from ..proj.reproj import reprojImg
from ..maths.fillnodata import replace_nans #inpainting function (ie fill nodata)
from ..maths.fillnodata import replace_nans, set_nans #inpainting function (ie fill nodata)
from ..utils import XY as xy
from ..checkdeps import HAS_GDAL, HAS_PIL, HAS_IMGIO
from .. import settings
Expand Down Expand Up @@ -459,6 +459,14 @@ def fillNodata(self):
# Inpainting
self.data = replace_nans(self.data, max_iter=5, tolerance=0.5, kernel_size=2, method='localmean')

def setNodata(self, setNodata=1.0):
# Cast to float
self.cast2float()
# Fill mask with NaN (warning NaN is a special value for float arrays only)
self.data = np.ma.filled(self.data, np.NaN)
self.data = set_nans(self.data, setNodata)


def reproj(self, crs1, crs2, out_ul=None, out_size=None, out_res=None, sqPx=False, resamplAlg='BL'):
ds1 = self.toGDAL()
if not self.isGeoref:
Expand Down
41 changes: 41 additions & 0 deletions core/maths/fillnodata.py
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,47 @@
DTYPEf = np.float32
#DTYPEi = np.int32

def set_nans(array, value):
"""
Replace NaN elements with a given value.

Parameters
----------
array : 2d np.ndarray
an array containing NaN elements that have to be replaced

value : float
Value with which to replace the NaN elements

Returns
-------
filled : 2d np.ndarray
a copy of the input array, where NaN elements have been replaced.
"""

filled = np.empty( [array.shape[0], array.shape[1]], dtype=DTYPEf)
print("set nans to %s" %value)

# indices where array is NaN
inans, jnans = np.nonzero( np.isnan(array) )

# number of NaN elements
n_nans = len(inans)

# fill new array with input elements
for i in range(array.shape[0]):
for j in range(array.shape[1]):
filled[i,j] = array[i,j]

for k in range(n_nans):
i = inans[k]
j = jnans[k]

# replace with value
filled[i,j] = value

return filled


def replace_nans(array, max_iter, tolerance, kernel_size=1, method='localmean'):
"""
Expand Down
18 changes: 16 additions & 2 deletions operators/io_import_georaster.py
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@
from ..core.proj import Reproj

from bpy_extras.io_utils import ImportHelper #helper class defines filename and invoke() function which calls the file selector
from bpy.props import StringProperty, BoolProperty, EnumProperty, IntProperty
from bpy.props import StringProperty, BoolProperty, EnumProperty, IntProperty, FloatProperty
from bpy.types import Operator

PKG, SUBPKG = __package__.split('.', maxsplit=1)
Expand Down Expand Up @@ -85,6 +85,11 @@ def listPredefCRS(self, context):
name="Specifiy raster CRS",
description="Specifiy raster CRS if it's different from scene CRS",
default=False )

setNodata: BoolProperty(
name="Set nodata values",
description="Assign value to existing nodata values to get an usuable displacement texture",
default=False )

# List of operator properties, the attributes will be assigned
# to the class instance from the operator settings before calling.
Expand Down Expand Up @@ -139,6 +144,12 @@ def listSubdivisionModes(self, context):
default=False
)
#
setNovalue: FloatProperty(
name="New nodata value",
description="Assign value to existing nodata values to get an usuable displacement texture",
default=0.0
)
#
step: IntProperty(name = "Step", default=1, description="Pixel step", min=1)

buildFaces: BoolProperty(name="Build faces", default=True, description='Build quad faces connecting pixel point cloud')
Expand Down Expand Up @@ -175,6 +186,9 @@ def draw(self, context):
if self.subdivision == 'mesh':
layout.prop(self, 'step')
layout.prop(self, 'fillNodata')
layout.prop(self, 'setNodata')
if self.setNodata:
layout.prop(self, 'setNovalue')
#
if self.importMode == 'DEM_RAW':
layout.prop(self, 'buildFaces')
Expand Down Expand Up @@ -384,7 +398,7 @@ def execute(self, context):

# Load raster
try:
grid = bpyGeoRaster(filePath, subBoxGeo=subBox, clip=self.clip, fillNodata=self.fillNodata, useGDAL=HAS_GDAL, raw=True)
grid = bpyGeoRaster(filePath, subBoxGeo=subBox, clip=self.clip, fillNodata=self.fillNodata, setNodata=self.setNodata, setNovalue=self.setNovalue, useGDAL=HAS_GDAL, raw=True)
except IOError as e:
log.error("Unable to open raster", exc_info=True)
self.report({'ERROR'}, "Unable to open raster, check logs for more infos")
Expand Down
50 changes: 46 additions & 4 deletions operators/mesh_earth_sphere.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
from bpy.types import Operator
from bpy.props import IntProperty

from math import cos, sin, radians, sqrt
from math import cos, sin, radians, sqrt, atan2, asin, degrees
from mathutils import Vector

import logging
Expand All @@ -13,9 +13,18 @@ def lonlat2xyz(R, lon, lat):
lon, lat = radians(lon), radians(lat)
x = R * cos(lat) * cos(lon)
y = R * cos(lat) * sin(lon)
z = R *sin(lat)
z = R * sin(lat)
return Vector((x, y, z))

def xyz2lonlat(x, y, z, r=100, z_scale=1):
R = sqrt(x**2 + y**2 + z**2)
lon_rad = atan2(y, x)
lat_rad = asin(z/R)
R = (R-r) / z_scale
lon = degrees(lon_rad)
lat = degrees(lat_rad)
return Vector((lon, lat, R))


class OBJECT_OT_earth_sphere(Operator):
bl_idname = "earth.sphere"
Expand All @@ -24,6 +33,7 @@ class OBJECT_OT_earth_sphere(Operator):
bl_options = {"REGISTER", "UNDO"}

radius: IntProperty(name = "Radius", default=100, description="Sphere radius", min=1)
z_scale: IntProperty(name = "Z scale", default=1, description="Scale for the z axis")

def execute(self, context):
scn = bpy.context.scene
Expand All @@ -50,8 +60,39 @@ def execute(self, context):
m = obj.matrix_world
for vertex in mesh.vertices:
co = m @ vertex.co
lon, lat = co.x, co.y
vertex.co = m.inverted() @ lonlat2xyz(self.radius, lon, lat)
lon, lat, z = co.x, co.y, co.z
vertex.co = m.inverted() @ lonlat2xyz(self.radius + z*self.z_scale, lon, lat)

return {'FINISHED'}

class OBJECT_OT_earth_plane(Operator):
bl_idname = "earth.plane"
bl_label = "sphere to latlon"
bl_description = "Transform data on a sphere to longitude/latitude on a plane"
bl_options = {"REGISTER", "UNDO"}

radius: IntProperty(name = "Radius", default=100, description="Sphere radius", min=1)
z_scale: IntProperty(name = "Z scale", default=1, description="Scale for the z axis")

def execute(self, context):
scn = bpy.context.scene
objs = bpy.context.selected_objects

if not objs:
self.report({'INFO'}, "No selected object")
return {'CANCELLED'}

for obj in objs:
if obj.type != 'MESH':
log.warning("Object {} is not a mesh".format(obj.name))
continue

mesh = obj.data
m = obj.matrix_world
for vertex in mesh.vertices:
co = m @ vertex.co
vertex.co = m.inverted() @ xyz2lonlat(co.x, co.y, co.z, self.radius, self.z_scale)


return {'FINISHED'}

Expand Down Expand Up @@ -92,6 +133,7 @@ def execute(self, context):

classes = [
OBJECT_OT_earth_sphere,
OBJECT_OT_earth_plane,
OBJECT_OT_earth_curvature
]

Expand Down
6 changes: 5 additions & 1 deletion operators/utils/georaster_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -217,7 +217,7 @@ def setDisplacer(obj, rast, uvTxtLayer, mid=0, interpolation=False):

class bpyGeoRaster(GeoRaster):

def __init__(self, path, subBoxGeo=None, useGDAL=False, clip=False, fillNodata=False, raw=False):
def __init__(self, path, subBoxGeo=None, useGDAL=False, clip=False, fillNodata=False, setNodata=False, setNovalue=0.0, raw=False):

#First init parent class
GeoRaster.__init__(self, path, subBoxGeo=subBoxGeo, useGDAL=useGDAL)
Expand All @@ -230,6 +230,7 @@ def __init__(self, path, subBoxGeo=None, useGDAL=False, clip=False, fillNodata=F
if self.format not in ['GTiff', 'TIFF', 'BMP', 'PNG', 'JPEG', 'JPEG2000'] \
or (clip and self.subBoxGeo is not None) \
or fillNodata \
or setNodata \
or self.ddtype == 'int16':

#Open the raster as numpy array (read only a subset if we want to clip it)
Expand All @@ -245,6 +246,9 @@ def __init__(self, path, subBoxGeo=None, useGDAL=False, clip=False, fillNodata=F
#replace nodata with interpolated values
if fillNodata:
img.fillNodata()

if setNodata:
img.setNodata(setNovalue)

#save to a new tiff file on disk
filepath = os.path.splitext(self.path)[0] + '_bgis.tif'
Expand Down