@@ -84,6 +84,13 @@ def make_single_type_map(
8484 filtered_map = habitat_map == habitat_value
8585 filtered_map .to_geotiff (filter_map_path , parallelism = max_threads )
8686
87+ if pixel_scale is not None and target_projection is not None :
88+ reprojected_area = filtered_map .area .reproject (
89+ yg .MapProjection (target_projection , pixel_scale , - pixel_scale )
90+ )
91+ else :
92+ reprojected_area = habitat_map .area
93+
8794 with VsimemFile (f"/vsimem/warped_{ habitat_value } .tif" ) as warped_map_path :
8895 logger .info ("Projecting %s..." , habitat_value )
8996 gdal .Warp (
@@ -97,6 +104,8 @@ def make_single_type_map(
97104 xRes = pixel_scale ,
98105 yRes = ((0.0 - pixel_scale ) if pixel_scale is not None else pixel_scale ),
99106 resampleAlg = "average" ,
107+ outputBounds = (reprojected_area .left , reprojected_area .bottom ,
108+ reprojected_area .right , reprojected_area .top ),
100109 targetAlignedPixels = pixel_scale is not None ,
101110 warpOptions = [f'NUM_THREADS={ max_threads } ' ],
102111 warpMemoryLimit = available_mem ,
@@ -107,6 +116,21 @@ def make_single_type_map(
107116 logger .info ("Saving %s..." , habitat_value )
108117 filename = f"lcc_{ habitat_value } .tif"
109118 with yg .read_raster (warped_map_path ) as result :
119+ # We had issues whereby original GDAL warp worked okay without needing
120+ # output bounds specified, and then at some point that changed and we were
121+ # losing a lot of the top and bottom of the map. This is a sanity check
122+ # just to ensure that isn't happening.
123+ # The 1% check here is probably overly tolerant, but it's enough to catch the
124+ # errors that caused this check to be added.
125+ reverted = result .area .reproject (habitat_map .map_projection )
126+ original = habitat_map .area
127+
128+ if ((abs (original .left - reverted .left ) / original .left ) > 0.01 ) or \
129+ ((abs (original .right - reverted .right ) / original .right ) > 0.01 ) or \
130+ ((abs (original .top - reverted .top ) / original .top ) > 0.01 ) or \
131+ ((abs (original .bottom - reverted .bottom ) / original .bottom ) > 0.01 ):
132+ raise ValueError (f"Area of reprojected map significantly different: { original } vs { reverted } " )
133+
110134 result .to_geotiff (output_directory_path / filename )
111135
112136def habitat_process (
@@ -145,7 +169,7 @@ def main() -> None:
145169 parser .add_argument (
146170 "--scale" ,
147171 type = float ,
148- required = True ,
172+ required = False ,
149173 dest = "pixel_scale" ,
150174 help = "Optional output pixel scale value, otherwise same as source."
151175 )
@@ -174,6 +198,11 @@ def main() -> None:
174198 )
175199 args = parser .parse_args ()
176200
201+ has_projection = args .target_projection is not None
202+ has_scale = args .pixel_scale is not None
203+ if has_projection ^ has_scale :
204+ parser .error ("Specify either both --projection and --scale must be specified or neither" )
205+
177206 habitat_process (
178207 args .habitat_path ,
179208 args .pixel_scale ,
0 commit comments