1. 通过GDAL warp来进行剪裁
import gdal
import os
in_shape = r'D:\shape_cn\guo.shp'
in_dir = r'D:\TEST'
out_dir = r'D:\TEST\output'
file_list = os.listdir(in_dir)
for file in file_list:
if file.endswith('.tif'):
print('Processing >>> ' + file)
in_raster = gdal.Open(os.path.join(in_dir, file))
out_raster = os.path.join(out_dir, file)
ds = gdal.Warp(out_raster,
in_raster,
format = 'GTiff',
cutlineDSName = in_shape,
# cutlineWhere="FIELD = 'whatever'", # clip specific feature
dstNodata = 0) # set nodata value
ds = None # close dataset
else:
print("No '.tif' file found ...")
print('..........ALL DONE..........')
2. 用rasterio来剪裁
import fiona
import rasterio
from rasterio.mask import mask
in_file = r'G:\test\test.tif'
shp = r"G:\Shapefiles\shape_cn\guo.shp"
with fiona.open(shp, "r", encoding='utf-8') as shapefile:
# obtain all the points of the polygon
geoms = [feature["geometry"] for feature in shapefile]
with rasterio.open(in_file) as src:
out_image, out_transform = mask(src, geoms, crop=True)
out_meta = src.meta.copy()
# update metadata for the output file
out_meta.update({"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform})
with rasterio.open(r"G:\test\masked.tif", "w", **out_meta) as dest:
dest.write(out_image)
print('... ... ... COMPLETED ... ... ...')
Python读取shapefile有以下几种方法:
- Fiona
import fiona
shape = fiona.open("my_shapefile.shp")
print shape.schema
{'geometry': 'LineString', 'properties': OrderedDict([(u'FID', 'float:11')])}
#first feature of the shapefile
first = shape.next()
print first # (GeoJSON format)
{'geometry': {'type': 'LineString', 'coordinates': [(0.0, 0.0), (25.0, 10.0), (50.0, 50.0)]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'FID', 0.0)])}
- PyShp
import shapefile
shape = shapefile.Reader("my_shapefile.shp")
#first feature of the shapefile
feature = shape.shapeRecords()[0]
first = feature.shape.__geo_interface__
print first # (GeoJSON format)
{'type': 'LineString', 'coordinates': ((0.0, 0.0), (25.0, 10.0), (50.0, 50.0))}
- ogr
from osgeo import ogr
file = ogr.Open("my_shapefile.shp")
shape = file.GetLayer(0)
#first feature of the shapefile
feature = shape.GetFeature(0)
first = feature.ExportToJson()
print first # (GeoJSON format)
{"geometry": {"type": "LineString", "coordinates": [[0.0, 0.0], [25.0, 10.0], [50.0, 50.0]]}, "type": "Feature", "properties": {"FID": 0.0}, "id": 0}
- geopandas
import geopandas as gpd
shapefile = gpd.read_file("shapefile.shp")
print(shapefile)