def __getitem__(self, geometry):
if isinstance(geometry, BaseGeometry) or getattr(geometry, "__geo_interface__", None) is not None:
if self._tms_meta._bounds is None:
return self.aoi(geojson=mapping(geometry), from_proj=self.proj)
image = GeoImage.__getitem__(self, geometry)
image._tms_meta = self._tms_meta
return image
else:
result = super(TmsImage, self).__getitem__(geometry)
image = super(TmsImage, self.__class__).__new__(self.__class__,
result.dask, result.name, result.chunks,
result.dtype, result.shape)
if all([isinstance(e, slice) for e in geometry]) and len(geometry) == len(self.shape):
xmin, ymin, xmax, ymax = geometry[2].start, geometry[1].start, geometry[2].stop, geometry[1].stop
xmin = 0 if xmin is None else xmin
ymin = 0 if ymin is None else ymin
xmax = self.shape[2] if xmax is None else xmax
ymax = self.shape[1] if ymax is None else ymax
g = ops.transform(self.__geo_transform__.fwd, box(xmin, ymin, xmax, ymax))
image.__geo_interface__ = mapping(g)
image.__geo_transform__ = self.__geo_transform__ + (xmin, ymin)
else:
image.__geo_interface__ = self.__geo_interface__
image.__geo_transform__ = self.__geo_transform__
image._tms_meta = self._tms_meta
return image
python类transform()的实例源码
def affine(self):
""" The geo transform of the image
Returns:
affine (dict): The image's affine transform
"""
# TODO add check for Ratpoly or whatevs
return self.__geo_transform__._affine
def _transpix(self, geometry, gsd, dem, proj):
xmin, ymin, xmax, ymax = geometry.bounds
x = np.linspace(xmin, xmax, num=int((xmax-xmin)/gsd))
y = np.linspace(ymax, ymin, num=int((ymax-ymin)/gsd))
xv, yv = np.meshgrid(x, y, indexing='xy')
if self.proj is None:
from_proj = "EPSG:4326"
else:
from_proj = self.proj
itfm = partial(pyproj.transform, pyproj.Proj(init=proj), pyproj.Proj(init=from_proj))
xv, yv = itfm(xv, yv) # if that works
if isinstance(dem, GeoImage):
g = box(xv.min(), yv.min(), xv.max(), yv.max())
try:
dem = dem[g].compute(get=dask.get) # read(quiet=True)
except AssertionError:
dem = 0 # guessing this is indexing by a 0 width geometry.
if isinstance(dem, np.ndarray):
dem = tf.resize(np.squeeze(dem), xv.shape, preserve_range=True, order=1, mode="edge")
return self.__geo_transform__.rev(xv, yv, z=dem, _type=np.float32)[::-1]
def _reproject(self, geometry, from_proj=None, to_proj=None):
if from_proj is None:
from_proj = self._default_proj
if to_proj is None:
to_proj = self.proj if self.proj is not None else "EPSG:4326"
tfm = partial(pyproj.transform, pyproj.Proj(init=from_proj), pyproj.Proj(init=to_proj))
return ops.transform(tfm, geometry)
def __contains__(self, g):
geometry = ops.transform(self.__geo_transform__.rev, g)
img_bounds = box(0, 0, *self.shape[2:0:-1])
return img_bounds.contains(geometry)
def get_tile_geometry(path, origin_espg, tolerance=500):
""" Calculate the data and tile geometry for sentinel-2 tiles """
with rasterio.open(path) as src:
# Get tile geometry
b = src.bounds
tile_shape = Polygon([(b[0], b[1]), (b[2], b[1]), (b[2], b[3]), (b[0], b[3]), (b[0], b[1])])
tile_geojson = mapping(tile_shape)
# read first band of the image
image = src.read(1)
# create a mask of zero values
mask = image == 0.
# generate shapes of the mask
novalue_shape = shapes(image, mask=mask, transform=src.affine)
# generate polygons using shapely
novalue_shape = [Polygon(s['coordinates'][0]) for (s, v) in novalue_shape]
if novalue_shape:
# Make sure polygons are united
# also simplify the resulting polygon
union = cascaded_union(novalue_shape)
# generates a geojson
data_shape = tile_shape.difference(union)
# If there are multipolygons, select the largest one
if data_shape.geom_type == 'MultiPolygon':
areas = {p.area: i for i, p in enumerate(data_shape)}
largest = max(areas.keys())
data_shape = data_shape[areas[largest]]
# if the polygon has interior rings, remove them
if list(data_shape.interiors):
data_shape = Polygon(data_shape.exterior.coords)
data_shape = data_shape.simplify(tolerance, preserve_topology=False)
data_geojson = mapping(data_shape)
else:
data_geojson = tile_geojson
# convert cooridnates to degrees
return (to_latlon(tile_geojson, origin_espg), to_latlon(data_geojson, origin_espg))
def getNoSampleGrid(yespoints, xvar, yvar, dx, h1, h2):
"""Return the grid from which "no" pixels can successfully be sampled.
Args:
yespoints: Sequence of (x,y) points (meters) where
landslide/liquefaction was observed.
xvar: Numpy array of centers of columns of sampling grid.
yvar: Numpy array of centers of rows of sampling grid.
dx: Sampling resolution in x and y (meters).
h1: Minimum buffer size for sampling non-hazard points.
h2: Maximum buffer size for sampling non-hazard points.
Returns:
Grid of shape (len(yvar),len(xvar)) where 1's represent pixels from
which "no" values can be sampled.
"""
shp = (len(xvar), len(yvar))
west = xvar.min() - dx / 2.0 # ??
north = yvar.max() + dx / 2.0 # ??
affine = affine_from_corner(west, north, dx, dx)
donuts = []
holes = []
for h, k in yespoints:
donut = createCirclePolygon(h, k, h2, dx)
hole = createCirclePolygon(h, k, h1, dx)
donuts.append(donut)
holes.append(hole)
donutburn = ((mapping(g), 1) for g in donuts)
holeburn = ((mapping(g), 2) for g in holes)
# we only want those pixels set where the polygon encloses the center point
alltouched = False
donutimg = rasterio.features.rasterize(donutburn,
out_shape=shp,
transform=affine,
all_touched=alltouched)
holeimg = rasterio.features.rasterize(holeburn,
out_shape=shp,
transform=affine,
all_touched=alltouched)
holeimg[holeimg == 0] = 1
holeimg[holeimg == 2] = 0
sampleimg = np.bitwise_and(donutimg, holeimg)
return sampleimg
def plot_predictions_for_a_city(df, name_of_predictions_col, city):
'''
INPUT: DataFrame with location predictions; name of column in DataFrame that
contains the predictions; city ('City, State') to plot predictions for
OUTPUT: Bokeh map that shows the actual location of all the users predicted to
be in the selected city
'''
df_ = df[df[name_of_predictions_col] == city]
# Initialize two lists to hold all the latitudes and longitudes
all_lats = []
all_longs = []
# Pull all latitudes in 'centroid' column append to all_lats
for i in df_['centroid']:
all_lats.append(i[0])
# Pull all longitudes in 'centroid' column append to all_longs
for i in df_['centroid']:
all_longs.append(i[1])
# Initialize two lists to hold all the latitudes and longitudes
# converted to web mercator
all_x = []
all_y = []
# Convert latittudes and longitudes to web mercator x and y format
for i in xrange(len(all_lats)):
pnt = transform(
partial(
pyproj.transform,
pyproj.Proj(init='EPSG:4326'),
pyproj.Proj(init='EPSG:3857')),
Point(all_longs[i], all_lats[i]))
all_x.append(pnt.x)
all_y.append(pnt.y)
p = base_plot()
p.add_tile(STAMEN_TERRAIN)
p.circle(x=all_x, y=all_y, line_color=None, fill_color='#380474', size=15, alpha=.5)
output_file("stamen_toner_plot.html")
show(p)
def trilaterate(points):
LatA = points[0][1]
LonA = points[0][0]
DistA = points[0][2]*MagicDistanceFactor
LatB = points[1][1]
LonB = points[1][0]
DistB = points[1][2]*MagicDistanceFactor
LatC = points[2][1]
LonC = points[2][0]
DistC = points[2][2]*MagicDistanceFactor
# Transform from Latitude/Longitude to ECEF coordinates.
xA,yA,zA = pyproj.transform(lla, ecef, LonA, LatA, 0, radians=False)
xB,yB,zB = pyproj.transform(lla, ecef, LonB, LatB, 0, radians=False)
xC,yC,zC = pyproj.transform(lla, ecef, LonC, LatC, 0, radians=False)
# Convert to numpy arrays.
P1 = numpy.array([xA/1000.0, yA/1000.0, zA/1000.0])
P2 = numpy.array([xB/1000.0, yB/1000.0, zB/1000.0])
P3 = numpy.array([xC/1000.0, yC/1000.0, zC/1000.0])
# Sphere intersection from Wikipedia + Stackoverflow.
ex = (P2 - P1)/(numpy.linalg.norm(P2 - P1))
i = numpy.dot(ex, P3 - P1)
ey = (P3 - P1 - i*ex)/(numpy.linalg.norm(P3 - P1 - i*ex))
ez = numpy.cross(ex,ey)
d = numpy.linalg.norm(P2 - P1)
j = numpy.dot(ey, P3 - P1)
x = (pow(DistA,2) - pow(DistB,2) + pow(d,2))/(2*d)
y = ((pow(DistA,2) - pow(DistC,2) + pow(i,2) + pow(j,2))/(2*j)) - ((i/j)*x)
# Handle errors.
if pow(DistA,2) - pow(x,2) - pow(y,2) < 0:
return []
z = numpy.sqrt(pow(DistA,2) - pow(x,2) - pow(y,2))
lon,lat,altitude = pyproj.transform(ecef, lla, x*1000,y*1000,z*1000, radians=False)
lon,lat,altitude = pyproj.transform(ecef, lla, x*1000,y*1000,z*1000, radians=True)
triPt = P1 + x*ex + y*ey + z*ez
lon,lat,altitude = pyproj.transform(ecef, lla, triPt[0]*1000,triPt[1]*1000,triPt[2]*1000, radians=False)
return [lat,lon]
# This was copied from trilat_optproblem.py and changed for only three point
# combinations.
def reproject_dataset_example(dataset, dataset_example, method=1):
# open dataset that must be transformed
try:
if dataset.split('.')[-1] == 'tif':
g = gdal.Open(dataset)
else:
g = dataset
except:
g = dataset
epsg_from = Get_epsg(g)
# open dataset that is used for transforming the dataset
try:
if dataset_example.split('.')[-1] == 'tif':
gland = gdal.Open(dataset_example)
else:
gland = dataset_example
except:
gland = dataset_example
epsg_to = Get_epsg(gland)
# Set the EPSG codes
osng = osr.SpatialReference()
osng.ImportFromEPSG(epsg_to)
wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(epsg_from)
# Get shape and geo transform from example
geo_land = gland.GetGeoTransform()
col=gland.RasterXSize
rows=gland.RasterYSize
# Create new raster
mem_drv = gdal.GetDriverByName('MEM')
dest1 = mem_drv.Create('', col, rows, 1, gdal.GDT_Float32)
dest1.SetGeoTransform(geo_land)
dest1.SetProjection(osng.ExportToWkt())
# Perform the projection/resampling
if method is 1:
gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_NearestNeighbour)
if method is 2:
gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Bilinear)
if method is 3:
gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Lanczos)
if method is 4:
gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Average)
return(dest1)
def __iter__(self):
""" Returns generator over shapefile rows.
Note:
The first column is an id field, taken from the id value of each shape
The middle values are taken from the property_schema
The last column is a string named geometry, which has the wkt value, the type is geometry_type.
"""
# These imports are nere, not at the module level, so the geo
# support can be an extra
self.start()
vfs, shp_file, layer_index = self._open_file_params()
with fiona.open(shp_file, vfs=vfs, layer=layer_index) as source:
if source.crs.get('init') != 'epsg:4326':
# Project back to WGS84
project = partial(pyproj.transform,
pyproj.Proj(source.crs, preserve_units=True),
pyproj.Proj(from_epsg('4326'))
)
else:
project = None
yield self.headers
for i,s in enumerate(source):
row_data = s['properties']
shp = asShape(s['geometry'])
row = [int(s['id'])]
for col_name, elem in row_data.items():
row.append(elem)
if project:
row.append(transform(project, shp))
else:
row.append(shp)
yield row
self.finish()
def read_records(records_csv, road_projection, record_projection, tz, col_id,
col_x, col_y, col_occurred):
"""Reads records csv, projects points, and localizes datetimes
:param records_csv: Path to CSV containing record data
:param road_projection: Projection CRS for road data
:param record_projection: Projection CRS for record data
:param tz: Timezone id for record data
:param col_id: Record id column name
:param col_x: Record x-coordinate column name
:param col_y: Record y-coordinate column name
:param col_occurred: Record occurred datetime column name
"""
# Create a function for projecting a point
project = partial(
pyproj.transform,
pyproj.Proj(record_projection),
pyproj.Proj(road_projection)
)
records = []
min_occurred = None
max_occurred = None
with open(records_csv, 'rb') as records_file:
csv_reader = csv.DictReader(records_file)
for row in csv_reader:
# Collect min and max occurred datetimes, as they'll be used later on
try:
parsed_dt = parser.parse(row[col_occurred])
# Localize datetimes that aren't timezone-aware
occurred = parsed_dt if parsed_dt.tzinfo else tz.localize(parsed_dt)
except:
# Skip the record if it has an invalid datetime
continue
if not min_occurred or occurred < min_occurred:
min_occurred = occurred
if not max_occurred or occurred > max_occurred:
max_occurred = occurred
records.append({
'id': row[col_id],
'point': transform(project, Point(float(row[col_x]), float(row[col_y]))),
'occurred': occurred
})
return records, min_occurred, max_occurred
def getGeonym(self, req, resp, query=None):
resp.status = falcon.HTTP_200
geo = None
# projections utilisées pour transformation en WGS84/Lambert93
s_srs = Proj(init='EPSG:2154')
t_srs = Proj(init='EPSG:4326')
if 'x' in req.params and 'y' in req.params:
lon,lat = transform(s_srs,t_srs,req.params['x'],req.params['y'])
query = geonym.ll2geonym(lat, lon)
elif 'lat' in req.params and 'lon' in req.params:
query = geonym.ll2geonym(float(req.params['lat']), float(req.params['lon']))
elif 'geonym' in req.params:
query = req.params['geonym']
elif 'adresse' in req.params:
r = requests.get(URL_GEOCODER+'/search', params={"q":req.params['adresse'], "autocomplete":0, "limit":1}, timeout=1)
geo = json.loads(r.text)
geo['source']=URL_GEOCODER
query = geonym.ll2geonym(geo['features'][0]['geometry']['coordinates'][1], geo['features'][0]['geometry']['coordinates'][0])
if query is not None and geonym.checkGeonym(query):
rev = None
data = geonym.geonym2ll(query)
if 'reverse' in req.params and req.params['reverse']=='yes':
r = requests.get(URL_GEOCODER+'/reverse', params={"lat":data['lat'],"lon":data['lon'],"limit":1}, timeout=1)
if r.status_code == 200:
rev = json.loads(r.text)
rev['source']=URL_GEOCODER
x,y = transform(t_srs,s_srs,data['lon'],data['lat'])
# on ne retourne les coordonnées Lambert que si on est en zone Lambert93
if y > -357823.2365 and x > 6037008.6939 and y < 1313632.3628 and x< 7230727.3772:
data['x'] = int(x)
data['y'] = int(y)
data['checksum'] = geonym.checksum(query)
geojson = {"type":"Feature",
"properties":data,
"link": "http://www.geonym.fr/visu/?g=%s" % (geonym.cleanGeonym(query),),
"params":geonym.getParams(),
"geometry":{"type":"Polygon","coordinates":[[[data['west'],data['south']],[data['east'],data['south']],[data['east'],data['north']],[data['west'],data['north']],[data['west'],data['south']]]]}}
if rev is not None:
geojson['reverse'] = rev
if geo is not None:
geojson['geocode'] = geo
resp.body = json.dumps(geojson, sort_keys=True, indent=4, separators=(',', ': '))
resp.set_header('Content-type','application/json')
else:
geojson = {
"type": "Feature",
"link": "https://github.com/geonym/geonymapi",
"params": geonym.getParams()
}
resp.status = falcon.HTTP_400
resp.set_header('Content-type', 'application/json')
resp.body = json.dumps(
geojson, sort_keys=True, indent=4, separators=(',', ': '))
