11 from typing
import List
13 PROGRAM =
"georegion_gen"
15 running_time = datetime.datetime.now()
16 RUNNING_TIME = running_time.strftime(
"%Y-%m-%dT%H%M%S")
20 with open(path_wkt)
as f:
21 wkts = f.read().splitlines()
22 return geopandas.GeoSeries.from_wkt(wkts)
29 for inp_shapefile
in inp_shapefiles:
31 for geom
in shape_data.geometry:
34 for inp_shapefile
in inp_shapefiles:
35 shape_data = geopandas.read_file(inp_shapefile)
36 for geom
in shape_data.geometry:
38 poly = geopandas.GeoSeries(poly)
39 west_lon = np.min(poly.bounds[
"minx"][:])
40 east_lon = np.max(poly.bounds[
"maxx"][:])
41 north_lat = np.max(poly.bounds[
"maxy"][:])
42 south_lat = np.min(poly.bounds[
"miny"][:])
43 if "global_set" in kwargs:
44 if kwargs[
"global_set"]:
51 xsize =
int((east_lon - west_lon) / scale_x) + 1
52 ysize =
int((north_lat - south_lat) / scale_y) + 1
53 shift_x_array = np.linspace(west_lon, east_lon, xsize)
54 shift_y_array = np.linspace(north_lat, south_lat, ysize)
55 trs = rasterio.transform.from_origin(west_lon, north_lat, scale_x, scale_y)
56 img = rasterio.features.rasterize(
57 poly, transform=trs, out_shape=(ysize, xsize))
59 f
"Generated georegion with shape {img.shape} along (lat,lon). Saving to {out_netcdf}")
60 with nc.Dataset(out_netcdf,
"w")
as georegion:
61 georegion.version = VERSION
62 georegion.program = PROGRAM
63 georegion.date_created = RUNNING_TIME
64 if "location" in kwargs
and not kwargs[
"location"]
is None:
65 georegion.location = kwargs[
"location"]
66 georegion.createDimension(
'lat', ysize)
67 georegion.createDimension(
'lon', xsize)
68 georegion_var = georegion.createVariable(varname=
'georegion', datatype=np.int8, dimensions=(
'lat',
'lon'), zlib=
True,
69 chunksizes=(432, 864), fill_value=-128)
70 lat_var = georegion.createVariable(
71 varname=
'lat', datatype=np.float64, dimensions=
'lat', zlib=
True, chunksizes=(432,), fill_value=-32767)
72 lon_var = georegion.createVariable(
73 varname=
'lon', datatype=np.float64, dimensions=
'lon', zlib=
True, chunksizes=(864,), fill_value=-32767)
74 georegion_var[:] = img
75 lat_var[:] = shift_y_array
76 lon_var[:] = shift_x_array
78 lon_var.standard_name =
"longitude"
79 lon_var.long_name =
"longitude"
80 lon_var.units =
"degrees_east"
83 lat_var.standard_name =
"latitude"
84 lon_var.long_name =
"latitude"
85 lat_var.units =
"degrees_north"
88 georegion_var.long_name =
"geographical_binary_mask"
89 georegion_var.description =
"Binary mask for masking out a specific earth domain "
90 georegion_var.comment =
"0 = outside, 1 = inside the region"
91 georegion_var.valid_min = 0
92 georegion_var.valid_max = 1
94 georegion.Conventions =
"CF-1.8 ACDD-1.3"
95 georegion.license =
"https://science.nasa.gov/earth-science/earth-science-data/data-information-policy/"
96 georegion.naming_authority =
"gov.nasa.gsfc.sci.oceandata"
97 georegion.keywords_vocabulary =
"NASA Global Change Master Directory (GCMD) Science Keywords"
98 georegion.standard_name_vocabulary =
"CF Standard Name Table v79"
99 georegion.creator_name =
"NASA/GSFC/OBPG"
100 georegion.creator_email =
"data@oceancolor.gsfc.nasa.gov"
101 georegion.creator_url =
"https://oceancolor.gsfc.nasa.gov"
102 georegion.project =
"Ocean Biology Processing Group"
103 georegion.publisher_name =
"NASA/GSFC/OB.DAAC"
104 georegion.publisher_email =
"data@oceancolor.gsfc.nasa.gov"
105 georegion.publisher_url =
"https://oceancolor.gsfc.nasa.gov"
106 georegion.institution =
"NASA Goddard Space Flight Center, Ocean Biology Processing Group"
107 georegion.product_name = os.path.basename(out_netcdf)
108 georegion.history =
" ".join(sys.argv)
109 georegion.geospatial_lat_min = south_lat
110 georegion.geospatial_lat_max = north_lat
111 georegion.geospatial_lon_min = west_lon
112 georegion.geospatial_lon_max = east_lon
113 georegion.geospatial_lat_units =
"degrees_north"
114 georegion.geospatial_lon_units =
"degrees_east"
115 georegion.geospatial_lat_resolution = scale_y
116 georegion.geospatial_lon_resolution = scale_x
117 georegion.grid_mapping_name =
"latitude_longitude"
121 if __name__ ==
"__main__":
122 print(PROGRAM, VERSION)
123 parser = argparse.ArgumentParser()
125 "--ifiles", help=
"input shape files", type=str, required=
True, nargs=
'+')
127 "--ofile", help=
"output georegion file", type=str, required=
True)
129 "--stepsize", help=
"stepsize along latitude", type=float, default=1.0 / 16.0 / 16.0)
133 "--pversion", help=
"pversion", type=str, default=
"V11")
135 "--global_set", help=
"global coverage", type=bool, default=
False)
137 "--wkt", help=
"indicate that input files are ASCII files which contains WKT strings", type=bool, default=
False)
138 if len(sys.argv) == 1:
141 print(
"You are running {}, version {}. Current time is {}".format(
142 PROGRAM, VERSION, running_time))
143 args = parser.parse_args()
145 shapefile_path = args.ifiles
146 stepsize = args.stepsize
149 global_set=args.global_set, wkt=args.wkt)
