Commit 2c88efe8 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Bugfix within GeoArray; added functions

geo.coord_trafo:
- added functions lonlat_to_pixel(), reproject_shapelyPoly()

geo.projection:
- added function get_prjLonLat()

io.raster.GeoArray:
- GeoArray.get_mapPos(): bugfix for returning the wrong target pixel size
- get_array_at_mapPos():  bugfix for returning the wrong target pixel size
parent 903d411c
......@@ -2,6 +2,8 @@
__author__ = "Daniel Scheffler"
import pyproj
import numpy as np
from shapely.geometry import Polygon
try:
from osgeo import ogr
......@@ -17,7 +19,7 @@ except ImportError:
from gdalconst import *
from .projection import get_proj4info
from .projection import get_proj4info, proj4_to_dict
......@@ -45,8 +47,8 @@ def transform_wgs84_to_utm(lon, lat, zone=None):
:param lat:
:param zone:
"""
assert -180.<=lon<=180. and -90.<=lat<=90., '''GEOPROCESSING error: Input coordinates for transform_wgs84_to_utm()
out of range. Got %s/%s.''' %(lon,lat)
assert -180.<=lon<=180. and -90.<=lat<=90., 'Input coordinates for transform_wgs84_to_utm() out of range. ' \
'Got %s/%s.' %(lon,lat)
get_utm_zone = lambda longitude: int(1 + (longitude + 180.0) / 6.0)
# def is_northern(latitude):
# """
......@@ -207,6 +209,17 @@ def latLonToPixel(latLonPairs, path_im=None, geotransform=None, projection=None)
return pixelPairs
def lonlat_to_pixel(lon, lat, inverse_geo_transform):
"""Translates the given lon, lat to the grid pixel coordinates in data array (zero start)"""
# transform to utm
utm_x, utm_y, utm_z = transform_wgs84_to_utm(lon, lat)
# apply inverse geo tranform
pixel_x, pixel_y = gdal.ApplyGeoTransform(inverse_geo_transform, utm_x, utm_y)
pixel_x = int(pixel_x) - 1 # adjust to 0 start for array
pixel_y = int(pixel_y) - 1 # adjust to 0 start for array
return pixel_x, abs(pixel_y) # y pixel is likly a negative value given geo_transform
def transform_GCPlist(gcpList, prj_src, prj_tgt):
"""
......@@ -217,4 +230,24 @@ def transform_GCPlist(gcpList, prj_src, prj_tgt):
:return:
"""
return [gdal.GCP(*transform_any_prj(prj_src, prj_tgt, GCP.GCPX, GCP.GCPY), GCP.GCPZ, GCP.GCPPixel, GCP.GCPLine)
for GCP in gcpList]
\ No newline at end of file
for GCP in gcpList]
def reproject_shapelyPoly(shapelyPoly, tgt_prj):
# type: (shapely.Polygon, str) -> shapely.Polygon
"""Repojects a shapely polygon that has LonLat coordinates into the given target projection.
:param shapelyPoly: <shapely.Polygon> the shapely polygon to be reprojected (must have lonlat coordinates)
:param tgt_prj: <str> WKT projection string (like GDAL projection)
"""
# TODO nicht sicher, ob wirklich nur lonlat funktioniert
get_coordsArr = lambda shpPoly: np.swapaxes(np.array(shpPoly.exterior.coords.xy), 0, 1)
coordsArr = get_coordsArr(shapelyPoly)
coordsArrTgtPrj = np.zeros_like(coordsArr)
proj4dict = proj4_to_dict(get_proj4info(proj=tgt_prj))
ReProj = pyproj.Proj(proj=proj4dict['proj'], zone=proj4dict['zone'], ellps=proj4dict['datum'])
x, y = ReProj(coordsArr[:, 0], coordsArr[:, 1], inverse=False)
coordsArrTgtPrj[:, 0] = x
coordsArrTgtPrj[:, 1] = y
return Polygon(coordsArrTgtPrj)
\ No newline at end of file
......@@ -120,4 +120,15 @@ def get_UTMzone(ds=None,prj=None):
srs.ImportFromWkt(ds.GetProjection() if ds else prj)
return srs.GetUTMZone()
else:
return None
\ No newline at end of file
return None
def get_prjLonLat(fmt='wkt'):
# type: (str) -> Any
"""Returns standard geographic projection (EPSG 4326) in the WKT or PROJ4 format.
:param fmt: <str> target format - 'WKT' or 'PROJ4'
"""
assert re.search('wkt',fmt,re.I) or re.search('Proj4',fmt,re.I), 'unsupported output format'
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
return srs.ExportToWkt() if re.search('wkt',fmt,re.I) else srs.ExportToProj4()
\ No newline at end of file
......@@ -378,7 +378,7 @@ class GeoArray(object):
'has to be passed.')
sub_arr, sub_gt, sub_prj = get_array_at_mapPos(self, arr_gt, arr_prj, mapBounds_prj, mapBounds, fillVal=fillVal,
rspAlg=rspAlg)
rspAlg=rspAlg, out_gsd=(self.xgsd,self.ygsd))
return sub_arr, sub_gt, sub_prj
......@@ -554,7 +554,7 @@ def get_array_at_mapPosOLD(arr, arr_gt, arr_prj, mapBounds, mapBounds_prj, band2
return out_arr, out_gt, out_prj
def get_array_at_mapPos(arr, arr_gt, arr_prj, out_prj, mapBounds, mapBounds_prj=None, band2get=None, fillVal=0,
def get_array_at_mapPos(arr, arr_gt, arr_prj, out_prj, mapBounds, mapBounds_prj=None, out_gsd=None, band2get=None, fillVal=0,
rspAlg='near'):
"""
......@@ -588,6 +588,6 @@ def get_array_at_mapPos(arr, arr_gt, arr_prj, out_prj, mapBounds, mapBounds_prj=
from ...geo.raster.reproject import warp_ndarray
out_arr, out_gt, out_prj = \
warp_ndarray(arr, arr_gt, arr_prj, out_prj=out_prj, out_bounds=mapBounds, out_bounds_prj=mapBounds_prj,
in_nodata=fillVal, out_nodata=fillVal, rspAlg=rspAlg)
in_nodata=fillVal, out_nodata=fillVal, rspAlg=rspAlg, out_gsd=out_gsd)
return out_arr, out_gt, out_prj
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment