fix white blank portions, fix georeference, fix NODATA
This commit is contained in:
parent
8fef56fa76
commit
a08a49611c
10
README.md
10
README.md
|
|
@ -1,3 +1,11 @@
|
|||
# wmts-scraper
|
||||
|
||||
Scrape tiles from WMTS servers.
|
||||
Scrape tiles from WMTS servers.
|
||||
|
||||
### ArcGIS
|
||||
|
||||
- If your raster is a teal square in ArcGIS, open Symbology and set `Stretch type` to `Esri`, which is the only stretch type that works with the background mask.
|
||||
|
||||
### Thank Yous
|
||||
|
||||
https://jimmyutterstrom.com/blog/2019/06/05/map-tiles-to-geotiff/
|
||||
|
|
@ -11,11 +11,9 @@ from rasterio import Affine
|
|||
from tqdm import tqdm
|
||||
|
||||
from pkg.image import random_file_width
|
||||
from pkg.spatial import deg2num
|
||||
from pkg.spatial import deg2num, lonlat_to_meters
|
||||
from pkg.thread import download_tile
|
||||
|
||||
# TODO: support reusing TIFFs
|
||||
|
||||
if __name__ == '__main__':
|
||||
parser = argparse.ArgumentParser(description='Exfiltrate data from WMS servers.')
|
||||
parser.add_argument('base_url', help='The base URL for the WMS server. Example: https://wmts.nlsc.gov.tw/wmts/nURBAN/default/EPSG:3857/')
|
||||
|
|
@ -82,11 +80,6 @@ if __name__ == '__main__':
|
|||
|
||||
def build_tiff_data(task):
|
||||
row, col = task
|
||||
# try:
|
||||
# row, col = map(int, file.name[:-4].split("_"))
|
||||
# except:
|
||||
# print('Bad file:', file)
|
||||
# return
|
||||
tile_file = tiles_output / f"{row}_{col}.png"
|
||||
if not tile_file.is_file():
|
||||
raise Exception(f'Tile does not exist: {tile_file}')
|
||||
|
|
@ -97,11 +90,20 @@ if __name__ == '__main__':
|
|||
# Remove the alpha channel
|
||||
tile_data = tile_data[:, :, :3]
|
||||
|
||||
# Calculate the position of the tile in the image data array
|
||||
# Replace white pixels with NODATA
|
||||
tile_data[np.all(tile_data == [255, 255, 255], axis=-1)] = [0, 0, 0]
|
||||
|
||||
# ArcGIS does not like pixels that have zeros in them, eg (255, 0, 0). We need to convert these to (255, 1, 1).
|
||||
mask = np.any(tile_data == 0, axis=-1) & np.any(tile_data != 0, axis=-1) # Identify pixels where not all bands are zero and at least one band is zero.
|
||||
for i in range(3): # Iterate over each band.
|
||||
# For these pixels, set zero bands to one.
|
||||
tile_data[mask & (tile_data[:, :, i] == 0), i] = 1
|
||||
|
||||
# Calculate the position of the tile in the image data array.
|
||||
row_pos = (row - min_row) * tile_size
|
||||
col_pos = (col - min_col) * tile_size
|
||||
|
||||
# Insert the tile data into the image data array
|
||||
# Insert the tile data into the image data array at the correct spot.
|
||||
image_data[row_pos:row_pos + tile_size, col_pos:col_pos + tile_size] = tile_data
|
||||
|
||||
|
||||
|
|
@ -110,18 +112,22 @@ if __name__ == '__main__':
|
|||
for future in tqdm(as_completed(futures), total=len(futures), desc='Building TIFF'):
|
||||
pass
|
||||
|
||||
# Transpose the image data array to the format (bands, rows, cols)
|
||||
# Transpose the image data array to the format (bands, rows, cols).
|
||||
image_data = np.transpose(image_data, (2, 0, 1))
|
||||
|
||||
# Define the transformation from pixel coordinates to CRS coordinates
|
||||
transform = Affine.translation(top_left_lon, top_left_lat) * Affine.scale((bottom_right_lon - top_left_lon) / (num_cols * tile_size), (bottom_right_lat - top_left_lat) / (num_rows * tile_size))
|
||||
# Convert geographic coordinates to Web Mercator coordinates. Not 100% sure this is nessesary.
|
||||
top_left_mx, top_left_my = lonlat_to_meters(top_left_lon, top_left_lat)
|
||||
bottom_right_mx, bottom_right_my = lonlat_to_meters(bottom_right_lon, bottom_right_lat)
|
||||
|
||||
# Define the CRS
|
||||
crs = rasterio.crs.CRS.from_string("EPSG:3857")
|
||||
# Define the transformation from pixel coordinates to geographic coordinates, which is an Affine transformation that
|
||||
# maps pixel coordinates in the image to geographic coordinates on the Earth's surface.
|
||||
transform = (Affine.translation(top_left_lon, top_left_lat) # Create a translation transformation that shifts the image and set the origin of the image to the top-left corner of the bounding box.
|
||||
# Create a scaling transformation that scales the image in the x and y directions to convert the pixel coordinates of the image to the geographic coordinates of the bounding box.
|
||||
* Affine.scale((bottom_right_lon - top_left_lon) / image_data.shape[2], (bottom_right_lat - top_left_lat) / image_data.shape[1]))
|
||||
|
||||
# Write the image data to a GeoTIFF file
|
||||
print('Saving to:', output_tiff)
|
||||
start = time.time()
|
||||
with rasterio.open(output_tiff, "w", driver="GTiff", height=num_rows * tile_size, width=num_cols * tile_size, count=3, dtype=str(image_data.dtype), crs=crs, transform=transform, compress="DEFLATE") as dst:
|
||||
with rasterio.open(output_tiff, "w", driver="GTiff", height=num_rows * tile_size, width=num_cols * tile_size, count=3, dtype=str(image_data.dtype), crs='EPSG:4326', transform=transform, compress="DEFLATE", nodata=0) as dst:
|
||||
dst.write(image_data, indexes=[1, 2, 3])
|
||||
print(f'Saved in {int(time.time() - start)} seconds.')
|
||||
|
|
|
|||
|
|
@ -1,9 +1,14 @@
|
|||
from rasterio import Affine
|
||||
|
||||
import math
|
||||
|
||||
|
||||
def deg2num(lat_deg, lon_deg, zoom):
|
||||
"""
|
||||
Convert a lat/long coordinate to x and y numbers.
|
||||
:param lat_deg:
|
||||
:param lon_deg:
|
||||
:param zoom:
|
||||
:return:
|
||||
"""
|
||||
lat_rad = math.radians(lat_deg)
|
||||
n = 2.0 ** zoom
|
||||
xtile = int((lon_deg + 180.0) / 360.0 * n)
|
||||
|
|
@ -11,25 +16,9 @@ def deg2num(lat_deg, lon_deg, zoom):
|
|||
return xtile, ytile
|
||||
|
||||
|
||||
def calculate_transform(args, min_col, max_col, min_row, max_row):
|
||||
# Constants for Web Mercator (EPSG:3857)
|
||||
tile_size_in_pixels = 256
|
||||
earth_circumference = 40075016.686 # in meters
|
||||
|
||||
# Calculate the size of a tile in meters at the given zoom level
|
||||
tile_size_in_meters = earth_circumference / (2 ** args.zoom)
|
||||
|
||||
# Calculate the geographic coordinates of the bounding box
|
||||
min_x = min_col * tile_size_in_meters
|
||||
max_x = (max_col + 1) * tile_size_in_meters
|
||||
min_y = (2 ** args.zoom - max_row - 1) * tile_size_in_meters
|
||||
max_y = (2 ** args.zoom - min_row) * tile_size_in_meters
|
||||
|
||||
# Calculate the pixel size in x and y directions
|
||||
pixel_size_x = tile_size_in_meters / tile_size_in_pixels
|
||||
pixel_size_y = tile_size_in_meters / tile_size_in_pixels
|
||||
|
||||
# Define the transformation
|
||||
transform = Affine.translation(min_x, min_y) * Affine.scale(pixel_size_x, -pixel_size_y)
|
||||
|
||||
return transform
|
||||
def lonlat_to_meters(lon, lat):
|
||||
origin_shift = 2 * math.pi * 6378137 / 2.0
|
||||
mx = lon * origin_shift / 180.0
|
||||
my = math.log(math.tan((90 + lat) * math.pi / 360.0)) / (math.pi / 180.0)
|
||||
my = my * origin_shift / 180.0
|
||||
return mx, my
|
||||
|
|
|
|||
Loading…
Reference in New Issue