texts: Reprojection of satellite images mosaic to CRS with latitude and longitudes

texts:
I downloaded MODIS product MOD13A1 for India. Then mosaiced the granules in GeoTIFF format using rasterio Python library. I have the appropriate shapefile to mask the mosaiced raster.

After mosaicing granules:

mosaic shape: (1, 12001, 14401)

crs: CRS({'proj': 'utm', 'zone': 35, 'ellps': 'GRS80', 'units': 'm', 'no_defs': True}) (found using mosaic.crs of rasterio library)

meta of mosaic file:

Shapefile I am using (converted to png to upload here):

CRS of shapefile: [epsg:4326]

I want to mask the mosaic using shapefile above.
To bring them in same CRS I converted mosaic to EPSG:4326

Code to reproject mosaic:

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

dst_crs = 'EPSG:4326'

with rasterio.open('mark_2.tif') as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
    })

with rasterio.open('reproject_mark_2.tif', 'w', **kwargs) as dst:
    for i in range(1, src.count + 1):
        reproject(
            source=rasterio.band(src, i),
            destination=rasterio.band(dst, i),
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=transform,
            dst_crs=dst_crs,
            resampling=Resampling.nearest)

Code to mask mosaic using shapefile:

import fiona
import rasterio
import rasterio.mask

with fiona.open("india_shape_file/IND_adm0.shp", "r") as shapefile:
    features = [feature["geometry"] for feature in shapefile]

with rasterio.open('mosaic.tif') as src:
out_image, out_transform = rasterio.mask.mask(src, features,
                                                    crop=True)
out_meta = src.meta.copy()

But something is going wrong.

Reprojected mosaic:

Meta of reprojected file:

After masking:

Mosaicing Code:

import rasterio
from rasterio.merge import merge
from rasterio.plot import show
import glob
import os

dirpath = '37758/*'
search_criteria = 'MOD13A1_A2010001*.tif'
q = os.path.join(dirpath, search_criteria)

files = glob.glob(q)

mosaic_files = []

for fp in files:
    src = rasterio.open(fp)
    mosaic_files.append(src)

mosaic, out_trans = merge(mosaic_files)
out_meta = src.meta.copy()
out_meta.update({"driver": "GTiff",
                 "height": mosaic.shape[1],
                 "width": mosaic.shape[2],
                 "transform": out_trans,
                 "crs": "+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs
"})

with rasterio.open('mark_2.tif', "w", **out_meta) as dest:
    dest.write(mosaic)

Reprojection is wrong. Can anyone tell me what am I doing wrong ? I have tried the same in QGIS but result is same. Finally, I have to do this programmatically.


texts: Set Numeric Format for Class Breaks Labels?

texts:
Why doesn't this code change the numeric format of the labels generated when I build this renderer? When I execute the set_Break() this nnumeric format is correctly set in the renderer but I still get numbers that look like "110.400000" in the generated label.

INumericFormat numericFormat = new NumericFormatClass();
numericFormat.AlignmentOption = esriNumericAlignmentEnum.esriAlignLeft;
numericFormat.AlignmentWidth = 0;
numericFormat.RoundingOption = esriRoundingOptionEnum.esriRoundNumberOfDecimals;
numericFormat.RoundingValue = 2;
numericFormat.ShowPlusSign = true;
numericFormat.UseSeparator = true;
numericFormat.ZeroPad = false;
IClassBreaksUIProperties classBreaksUIProperties =
classBreaksRenderer as IClassBreaksUIProperties;
classBreaksUIProperties.NumberFormat = numericFormat as INumberFormat;
classBreaksRenderer = classBreaksUIProperties as IClassBreaksRenderer;


texts:
I'm trying to generate rasters of slope and curvature of an elevation DTM. The particularity is that i need to be able to specify the direction (e.g 180, 270 degrees)in this calculation.

I am using ArcGIS Desktop 10, however if the solution comes from another software package or code it's fine.


texts:
I just starting to work with raw-OSM data, so maybe my question is a very basic one.

When I download OSM data from Geofabrik I never get the same coverage of the polygons as in the rendered raster OSM data - at least for many areas.

In the screenshot below you see the rendered OSM data as background and the polygons of the downloaded raw-OSM data above. For some areas I get the polygons (the light-red colored area, transparent to see the rendered OSM data below), but you see e.g. the forest in the middle of the sceenshot and the farmland in the lower part which are missing in the raw-OSM data. Even the boundary between forest and agricultural area is not available (e.g. as Polyline in the .osm.pbf file).

Being used to work with shapefiles, I would assume that all areal features like land use are provided as polygons. But this is not the case obviously. The alternative -which I do not beliefe to be true - would be that the raw data have gaps in it.

Can someone explain what I am doing wrong?

Thanks

some technical details: I tried different data sets (meaning: different areas) and data formats (shapefiles, osm.pbf) but always have the same issue. I also tried to download OSM data via the OpenStreetMap plugin in QGIS.


texts: From gpx to line segment 1Km using ArcGIS 10

texts:
I have gpx files containg points with attribute data such as time and date. I'd like to convert them into lines, and that is not a problem, but what I want to do is create 1Km segments with time and date attribute.
I'm using Arcgis 10.2

Thanks


texts: QGIS crashes when loading sid file, gdalinfo also produces error in terminal

texts:
I have the same problem as this person QGIS crashes when I load .sid layer and I am running qgis v2.18 on mac os 10.12.4. I have gdal 1.11 and 2.1 installed and have also made sure to install the correct files in the libraries folders for mrsid (downloaded from lizard tech) per instructions on this site  http://www.kyngchaos.com/software/frameworks. Every time I try to add a sid file to qgis it crashes. when I try to run gdalinfo on the image in terminal I get this error:

Segmentation fault: 11

Does anyone know what might be wrong?


texts:
I am completely lost in the dozen of solution I've found. I would like an easy solution to create an OSM basemap for my different projects with nice styling.

I work directly on my computer, no server or web-based service.

I have set up a postgis database and imported my country's OSM extract (from geofrabick) in it.

And now I would like to know how can I use this in QGis as a basemap with mapnik style or other nice style? I don't need this basemap to be editable nor vertor. Idealy juste a tile raster map.

I know there is somme style sheet around, but if you use large OSM extract rendering if ver very long and not usable to work with you project.


texts: Can't run processing.algrun() in external Python application

texts:
I am trying to create an external Python application. I want to run this as a stand alone program in command line. I am trying to not use the QGIS console.

A first step is to clip a point shapefile using a polygon shapefile.

When I try to use processing.runalg("saga:clippointswithpolygons",...) nothing happens. There are no error messages or warnings. Nothing is written to the output shapefile.

import sys
import urllib

import pathlib

from qgis.core import QgsApplication
from qgis.core import QgsVectorLayer
from qgis.core import QgsVectorFileWriter

# Constants
node_shape_path = r'C:\Temp\TMC_nodes.shp'
clip_path = r'C:\Temp\clip.shp'  # Path to clipping polygon
out_shapefile_path = r"C:\Temp\TMC_Nodes_Clipped_2.shp"

##
# Prepare the QGIS environment and framework
##

# Prepare the QGIS environment
qgis_prefix = u'C:/PROGRA~1/QGISPI~1/apps/qgis'
app = QgsApplication([], True)
QgsApplication.setPrefixPath(qgis_prefix, True)
QgsApplication.initQgis()

# Prepare the QGIS processing framework
sys.path.append(plugins_path)
#from processing.core.Processing import Processing
import processing

##
# Clip the nodes and write to a new shapefile
##
processing.runalg("saga:clippointswithpolygons", node_shape_path, clip_path, "OBJECTID", 0, out_shapefile_path)

I have tried different versions of the processing.runalg line where I use file path strings instead of layers. I even went so far as to take the point layer, made from the csv, write it to a shapefile, and then use that path.

I have used the processing toolbox GUI to successfully conduct the clip. The log file reads:

ALGORITHM|Mon Sep 21 2015 12:26:24|
processing.runalg("saga:clippointswithpolygons","C:/Temp/TMC_nodes.shp","C:/Temp/FCMS_clip.shp","OBJECTID",0,"C:/Temp/TMC_Nodes_Clipped.shp")

I have even tried pasting the command call from the log directly into my external script, but it still does not work.


texts: Soil interpolation in ArcGIS 3D Analyst?

texts:
I have borehole data with different soil classifications.

I created each soil type in to individual shapefile.

Does anybody have any ideas on how to display a 3D subsoil model in ArcScene?

Can you tell me the steps to interpolate and finally visualize it in 3d with layers?

This screenshot is of my data.


texts: ArcMap layout view missing page rulers and border?

texts:
Sometimes, but not always, when I open an existing .mxd, the rulers and border (and shadow) are missing in layout view. All elements are there otherwise. When this happens, I am unable to export to either .emf or .pdf. I can't discern any rhyme or reason.

These are data driven pages, in case that's an issue.


texts: PostGIS 2.0 "No function matches the given name and argument types"

texts:
After a lot of research, I figured out how to build PostGIS 2.0 from source on a Turnkey Linux  LAPP (Postgres 8.4) server. The main draw of PostGIS 2.0 for me is the ST_ConcaveHull function, aka "shrink wrap".

Now that I have a functioning PostGIS 2.0, I'm trying to implement these functions.

I'm working with the NWS Shapefiles for States and Forecast Zones. Most if not all of the geometries I'm dealing with are multipolygons.

I've figured out how to draw these polygons, and this is the query I use to request data in the GeoJSON format:

SELECT ST_AsGeoJSON(geom) as geom FROM states WHERE id = 5;

Unfortounately, because I'm working with State sized data, some of the requests are HUGE! The GeoJSON for New Jersey is 2mb! I know PostGIS has some great functions to reduce the complexity of the geometries, but I'm having trouble implementing them.

The functions I want to use are:

ST_Simplify
ST_Union
ST_ConcaveHull

Unfortunately, I'm getting some very unhelpful error messages when I run queries using these functions.

SELECT ST_AsGeoJson(ST_Simplify(geom)) as geom FROM states WHERE id = 5;

Error: SQLSTATE[42883]: Undefined function: 7 ERROR: function st_simplify(geometry)
does not exist LINE 1: SELECT ST_AsGeoJson(ST_Simplify(geom)) as geom FROM states W...
^ HINT: No function matches the given name and argument types. You might need to add
explicit type casts.

I went through a LOT of time and effort to make sure that I had PostGIS 2.0 instead of 1.5 for the explicit purpose of using things like ST_ConcaveHull and ST_Simplify, and now Postgre is telling me they don't exist!

I found another post on here that gave me a query listing all global functions available to Postgre:

SELECT  proname, proargnames, prosrc
FROM    pg_catalog.pg_namespace n
JOIN    pg_catalog.pg_proc p
ON      pronamespace = n.oid
WHERE   nspname = 'public'

Listed in the results are:

st_union {geom1,geom2} geomunion
st_union pgis_union_geometry_array
st_union aggregate_dummy
st_concavehull {param_geom,param_pctconvex,param_allow_holes}
st_simplify LWGEOM_simplify2d
st_simplifypreservetopology topologypreservesimplify

So I'm kind of stuck and confused at this point. The query results say the function doesn't exist, but the output of the above query says they DO exist.

Can anybody shed some light on this for me? What am I doing wrong?


texts: Extracting all vertices from line shapefile in a simple way? (no python scripting)

texts: I found that shpdump is a very good tool to extract all vertices from an ESRI polyline shapefile, but I need a single triplet (X,Y,Z) file (ie. shpdump provides all coordinates by parts). Is there a similar command-line tool that can do this? I found some python solutions, but I need something like ogr2ogr or shpdump. ogr2ogr is a wonderful tool but it seems that I cannot extract vertices from a polyline shapefile.


texts: old coordinate systems that were "made up" by surveyors decades ago

texts:
Using ArcGIS COGO tool for ROW / Easement digitizing from old plats (using their bearings and distance calls) that used what a veteran of the county says were 'made up" coordinate systems in Maryland. Several cities supposedly had their own "made up" coordinate systems.

See this plat for example

Does this veteran's story sound correct?
I have contacted several people at state archives (source of the old plats) but no one there has answers to where / what these old coordinate systems are.

Is it possible in ArcGIS COGO or QGIS to recreate these "made up" coordinate systems, do the bearings / distance digitizing, and then transform them for use in Maryland State Plane 1900?


texts: Avoid Intersections not honored when saving in QGIS?

texts:
I'm editing a polygon layer in QGIS that contains several small polygons within larger polygons.  To save time drawing the larger (encompassing) polygons, I'm using the "avoid intersections" option, along with "enable topological editing" and "enable snapping on intersection" options on.  When I complete drawing a large polygon, all polygons render properly (i.e., the smaller polygons within the larger polygon are avoided/preserved as individual features, with proper topology as expected).  However, once I save my edits, the larger polygon covers the smaller polygons.  Is this an issue anyone else has encountered?  If so, what is the solution?

I have tried exporting the shapefile (save as--> new name).  That does not work.

The attached photos illustrate more clearly what I am trying to describe.  The image above shows the smaller polygons that were drawn first.  The image below shows the polygon I have drawn on top of them (colored yellow), which covers them upon saving the edits.


texts: Determining the area of each polygon within a polygon boundary?

texts:
I have a series of maps representing a 1000m diameter around a point. These maps contain 40 different types of polygons that represent a land use (ex. Corn crop, Pasture).

I want to select all polygons that are within the 1000m boundary and find out each of their areas.

I've attempted to create unions and delete the exterior polygons, but I'm only limited to performing a union on 2 features which would make that task very tedious. Is there a faster and simpler way of performing this task?

My goal is to create an attbiute table, that will have a unique identifer for each polygon and then a column with the SHAPE_AREA that I can export to excel.

I created the maps using ArcGIS, I have access to all licenses and other GIS programs.
All polygons are shapefiles, including the boundary circles.


texts: data type conversion extracting user defined tags as columns using Openstreetmap data

texts:
we have an ,osm data in which we have used some tags that we defined ourselves like, age, roof_type, roof_slope . we want to convert this file to shapefile, in which the roof_type tag appears as a column name and its values  as attributes.

we have tried:

QGIS: install plugin, only creates column of 7 tags, others are grouped into one column called tags

OSM2GIS: online tools, selects area on map, has many tags, but not user defined ones

Geoconverter: online tools, input a supported file and it will export a supported file, but not user defined ones

OSMembrane: still learning, frontend for Osmosis


texts: Get info from WMS layer with OpenLayers mobile

texts:
I am trying to build a web application for iPad using OpenLayers Mobile. One of the features it should have is that you can tap on a tree displayed on a basemap and get information about this tree. These trees are provided through a WMS service in GeoServer and the info should be displayed in a popup.

In the standard version of OpenLayers I would use WMSGetFeatureInfo to achieve this, but this seems not to be working with the mobile version.

When using OpenLayers mobile, I get an error message saying:

TypeError: OpenLayers.Control.WMSGetFeatureInfo is not a constructor

Can you help me fix this?


texts: Which spatial bundling or grouping algorithm optimally groups elements by attribute?

texts:
As I delve deeper into GIS, I find myself trying to improve efficiencies of certain procedures and processes.
The one process I'm working on now is spatially scheduling construction projects. For example, consider the two construction schedules outlined in the images below. The first one is perhaps the defacto non spatially planned schedule; however, the second, is more optimal and due to bundling of contracts would be cheaper as well.
This is easy when the network has 30 links; however, the one I am working on has ~50,000 links.
Question: What should I read upon or learn to implement such algorithms? Further more, are there any application type paper that you can reference to help my understanding on the matter?

Which GIS algorithm can I use to move the schedule from above to something similar as in the bottom image?

Note: I'm not looking for a software tool to do the work for me, I'm looking to understand the process and perhaps implement it as a plugin in QGIS or ArcGIS


texts:
I am using GDAL in Python to create a tiff raster of data that is originally from a shapefile. I have an auxillary data channel, which I would like to use to color the pixels of the tif. I have made sure that my data (yield_grid) is the same size and shape as my tif raster will be, but so far I am unable to get anything besides an all black square out of my raster.

Does anyone have a suggestion for what I might be doing wrong in assigning my data to my raster band?

# Define pixel_size and NoData value of new raster
pixel_size = 0.00005
NoData_value = -99

# Filename of input OGR file
vector_fn = inShapeFile

# Filename of the raster Tiff that will be created
raster_fn = 'test.tif'

# Open the data source and read in the extent
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
source_srs = source_layer.GetSpatialRef()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# Create the destination data source
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)

# Get Driver and Create output target
target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res,1,gdal.GDT_Byte)

# Set Geo Transform
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))

# Write array into raster band, normalize by max, and scale to 1024
target_ds.GetRasterBand(1).WriteArray(np.transpose(yield_grid/np.max(yield_grid))*1024)

# Set Spatial reference and projection
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromEPSG(4326)
target_ds.SetProjection(outRasterSRS.ExportToWkt())

#Flush Cache for raster band
target_ds.GetRasterBand(1).FlushCache()

# Printing raster band to make sure that I'm not crazy
print target_ds.GetRasterBand(1).ReadAsArray()
[[  0   0   0 ...,  39   2   0][  0   0   0 ...,  80  37   0]
 [  0   0   0 ...,  85  32   0]
 ...,
 [  0  58 119 ...,   0   0   0]
 [  0   0  82 ...,   0   0   0]
 [  0   0   0 ...,   0   0   0]]

The output at the end of the code shows that I am in fact putting values into the raster band, but I'm not sure why the tiff isn't showing any colors.


texts:
I would like to search for the specific buildings (polyhouse in my case) on the map like the one shown here

Are there any free options available which can give me lat and lng of such shapes/features in specified area?


texts: Adding own items into QGIS 3.10 news feed?

texts: Is there a possibility of adding own items into the news feed of QGIS 3.10? As a GIS administrator it would be a great tool to spread GIS/QGIS related news to internal colleagues.


texts: How to create SQL statement to only return chainage values every 500m?

texts:
I have a list of chainage values every 100m (as a whole number) and I am trying to create an SQL statement to only return values every 500m. Any ideas?

Using ArcGIS 10.1


texts: find route distance & time taken to travel between two points using leaflet

texts:
I am able to find route between two points using cloud mode & leaflet using this.
CloudMade Routing with Leaflet: route is not displayed
But I also want to find distance & time taken to travel from first point to last point.


texts:
The below Python code opens "select vector layer..." dialogbox. Selecting one or more layers works nicely.

How can we specify the layers programmatically using Python and then skipping the dialogbox? The dialogbox becomes annoying very quickly and it is the same layers that is needed.

qgis.utils.iface.addVectorLayer('C:\tmp\enc\NN1234.00','tbd','ogr')

Using QGis 2.18.15


texts: GeoServer WFS OpenLayers EditingPanel Problem

texts:
I'm trying to use OpenLayers to edit features of vector data served from a local geoserver. So far, I'm able to pull WMS and WFS layers using OpenLayers from the same geoserver. But now, I'd like enable editing of features of the vector layers. The relevant code section is as shown below:

var ke_admin_vector = new OpenLayers.Layer.Vector("Admin Vector", {
           strategies: [new OpenLayers.Strategy.BBOX(), new OpenLayers.Strategy.Save()],
           protocol: new OpenLayers.Protocol.WFS({
               version: "1.1.0",
               url: "http://HOST_IP:port/geoserver/wfs",
               featureType: "ke_admin",
               featureNS: "pop_explorer",
               srsName: "EPSG:4326"
           })
       });

       var editpanel = new EditingPanel(ke_admin_vector)
       map.addControl(editpanel);

ke_admin is the name of the layer. pop_explorer is the workspace url. EPSG:4326 is the projection of the layer within geoserver.

When I try accessing this from a browser, firebug complains that EditingPanel is not defined. Any suggestions on what I may be doing wrong. Thanks in advance.

ANSWERED.


texts:
Using Arc Desktop 10.3.1 I have a datetime field that is in the format of MM/DD/YY HH:MM:SS   Some of the records have a only a date and some have both a time and a date.  I want to strip off the time part using .split(" ")[0] but I am having difficulty figuring out how to craft the query to only select the records that have a time part.  I have tried:

select * from !DateTime!.split(" ")[1] > 0

The logic being I want to select all records from DateTime where the trailing time is greater than 0, then after they are selected I will run the field calculator to strip the time.


texts: OpenLayers: How can I add multiple KML files to the same vector layer?

texts:
I'm trying to add some KML files to a vector layer using this code:

            var map = new OpenLayers.Map('map');
            map.addLayer(new OpenLayers.Layer.OSM());
            var kml = ['kml/file1.kml','kml/file2.kml' ... 'kml/fileN.kml'];

            function createFeatures() {
                var features = [];
                for (var i=0; i<kml.length; i++) {
                        features.push(new OpenLayers.Feature.Vector(
                            new OpenLayers.Protocol.HTTP({
                                url: 'kml/file'+i+'.kml',
                                format: new OpenLayers.Format.KML()
                            }
                        )));
                }
                return features;
            }

            var vectors1 = new OpenLayers.Layer.Vector("Vector Layer 1", {
                isBaseLayer: false,
                style: { strokeWidth: 3, strokeColor: "#FF0000", fillOpacity: 0.50, fillColor: "#FF0000" }
            });
            map.addLayer(vectors1);
            map.addControl(new OpenLayers.Control.LayerSwitcher());
            vectors1.addFeatures(createFeatures());

I get this error:

Uncaught TypeError: a.geometry.getBounds is not a function


texts: assigning census block demographics to postal data (zip/zip+4)

texts: Using Python and R, what would be the best strategy to tie Census Block-level data to zip/zip + 4? I understand we're not matching polygon boundaries (zip codes are not geographies), but somehow service bureaus/lettershops tie addresses to Census demographic data. Do they work off zip+4 centroids? How can this lookup be updated with new zip/zip+4 definitions (i.e. not be bound by the 2010 ZTCA). Thanks for any pointers/links/book references.


texts: NAD 1983 UTM ZONE 10N or WGS 1984 UTM ZONE 10N?

texts:
I calculate small acreages of areas in Sacramento, no larger than about 2 acres at a time. I need an accurate projection to measure by the T as I don't always have access to groundtruthing.

Thoughts?


texts: Rendering raster layer from DTED based on height color map

texts:
I have a raster layer that was create from DTED file (dt1 specifically), created by the following code:

var rasterLayer = new RasterLayer();
rasterLayer.CreateFromFilePath(dtedFilePath);

The dted file is treated like raster and is displayed as grayscale.

From the other hand, I have a color map file (.clr) that maps specific colors to specific elevation values.
The file is in the form

0 63 187 239
1 63 187 239
2 125 198 34
...
2806 220 220 220

I would like, somehow, to have that coloring (the one in the clr file) used by some kind o renderer to render the aforementioned raster layer.

Is there any way to do this?

Thanx

Edit: forgot to mention, we are talking about a windows form application written in C#, with a simple map control in the center. I'm using ArcGIS Engine 10


texts: DEM water body visualisation (QGIS): extract areas with a min amount of pixels with the same value (= water surface)

texts:
I want to create a background map of a DEM (1m spatial resolution) to show the study area.
The study area (Northern Germany) contains many water bodies. I have problems with the visualisation of these water bodies (using QGIS) as

the lakes surfaces do have different pixel values, so it's not possible to just set a minimum threshold to visualize the water bodies because then areas of the non-water terrain fall under this category (info: each lake has one pixel value) --> see picture below

and

the lake-shapefiles that I have (e.g. OSM or ATKIS DLM) are not precise enough (measuring instruments are positioned right at the shores or in a certain distance from the shores within and outside the water bodies --> DEM shows the most precise boundaries).

I am thinking about extracting a mask with the water bodies from the DEM, but I am not sure about the method.
In my understanding, I could buffer (approx 2m) the lake-shapefiles that I have (they do mainly display the existing water bodies a little bit too small), and then extract areas with a min amount of pixels with the same value (= water surface) - but I don't know how that could possibly work...
Does anybody have a tip how to extract the water bodies or how to color areas with a certain amount of pixels with the same value next to each other (= the water bodies)?


texts: Store and retrieve date of geoproccessing in feature class

texts:
I have geo-processing models and scripts that are run on existing feature classes. The models and scripts truncate(delete) the rows in feature classes, and then append new rows. The geo-processing models and scripts are run manually/as-needed by a handful of users. The database is an Oracle Esri.st_geometry geodatabase, and the application is ArcGIS Desktop 10.3.1.

Is there a way to store the date that the geo-processing model or script is run in the feature class?

I don't want to store the date in the table name, because the feature class name is referenced by other systems; I wouldn't want the name to be constantly changing. I also need to be able to retrieve the date via some sort of script or query.


texts: Neighborhood Analysis with ArcGIS Spatial Analyst to Calculate Percentage

texts:
I have a raster file where the cells (1 meter) are classified as 1 (forest) and 0 (non-forest). I would like to create another raster (1 meter) where each cell is given the percentage of cells within 60 meters that are = 1 (forest).

I've looked at the Block Statistics and Focal Statistics tools but they do not quite give me what I'm looking for. The statistics are based on the values of the cells (1,0) not the count of cells in the search radius to calculate a percentage.

I'd like to keep this as a raster analysis rather than converting to points.

I'm using ArcGIS 10.5.1, Advanced.


texts: What is the unit of measurement of the distance between two points using geopandas distance function?

texts:
I have created three buffers using geopandas in python and I have a few points scattered here and there, some of them are inside the buffers.But the problem is the following:

When I obtain distance between the buffer centers and the points, for those points which are inside the buffer the distance shown here is 0 and for others some value >0. But the problem may arise if we have all these randomly scattered points inside the buffers, then the distance from the center of the buffer to all these points is 0. But in that case you cannot sort the values distance wise (in meters) and find the minimum one?

My code here is thus:

    import pandas as pd
    import geopandas
    import numpy as np
    us_states= geopandas.read_file("C:/Users/guptac/Downloads/cb_2017_us_state_500k/cb_2017_us_state_500k.shp")
    print(us_states.geometry[0])

    from shapely.geometry import Point,Polygon,LineString
    import fiona
    import shapely
    from shapely.geometry import shape

    df = pd.DataFrame(
    {'City': ['Fargo', 'Orange', 'Jersey City'],
     'State': ['North dakota', 'California', 'New Jersey'],
     'Latitude': [46.877186 , 33.787914 , 40.728157 ],
     'Longitude': [-96.789803, -117.853104, -74.077644]})
    df['Coordinates']  = list(zip(df.Longitude, df.Latitude))
    df['Coordinates'] = df['Coordinates'].apply(Point)
    gdf = geopandas.GeoDataFrame(df, geometry='Coordinates')
    gdf_points=gdf
    print(gdf.head())

    world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))
    #nybb=  geopandas.read_file(geopandas.datasets.get_path('nybb'))

    # We restrict to USA.
    ax = world[world.name == 'United States'].plot(
    color='white', edgecolor='black',figsize=(20,10))
    us_states.plot(ax=ax,figsize=(20,10))
    gdf.plot(ax=ax, color='red')
    ax.set(xlim=(-140,-50),ylim=(0,60))
    us_states.crs
    gdf.crs= {'+init' :'epsg:4326'}
    print(gdf.crs['units'])
    us_name = world[world.name=='United States']
    gdf['Coordinates']=gdf.buffer(3)
    #gdf.head(5)
    gplot=geopandas.overlay(us_states,gdf,how='identity')
    gx=gplot.plot(edgecolor='k', alpha=0.5, cmap='tab10', ax=ax)
    gdf.plot(ax=ax)

    df1=pd.DataFrame(
            {'Points':['Point1','Point2','Point3','Point4','Point5','Point6'],
             'Latitude':[44.773491,30.287765,41.279906,45.234565,31.23345,38.778987],
             'Longitude':[-94.789323,-114.565309,-72.055234,-93.230988,-112.233456,-72.239989]})
    df1['Coordinates']  = list(zip(df1.Longitude, df1.Latitude))
    df1['Coordinates'] = df1['Coordinates'].apply(Point)
    gdf1 = geopandas.GeoDataFrame(df1, geometry='Coordinates')
    print(gdf1.head())

    gdf1.plot(ax=gx,color='red')
    gx.set(xlim=(-140,-50),ylim=(0,60))

    #a = Point(-93.230988, 45.234565)
    centroid_of_buffers=gdf['Coordinates'].centroid
    #dist_a=centroid_of_buffers.distance(a)
    #min_dist=centroid_of_buffers.distance(gdf1[gdf1.Coordinates=="POINT (-93.230988 45.234565)"])
    gdf_distance_cluster_Nevada=gdf1.distance(gdf['Coordinates'][0],)
    gdf_distance_cluster_California=gdf1.distance(gdf['Coordinates'][1])
    gdf_distance_cluster_New_Jersey=gdf1.distance(gdf['Coordinates'][2])

From here what I can guess is how do I obtain the euclidean distance between the buffer centers and the  randomly scattered points ?

You can find the us_states dataset from https://www.census.gov/geo/maps-data/data/cbf/cbf_state.html


texts: Opening file geodatabases created with ogr2ogr in ArcGIS Desktop?

texts:
I've been using the new FileGDB support in gdal-1.9.0 to create file geodatabases that can be opened in ArcGIS.  For example (I'm doing this all on one line, so I'm not certain this will run as indented here):

ogr2ogr -overwrite -f 'FileGDB' \
    -sql "SELECT
              CAST(gml_id AS character(255)), id,
              CAST(scene_id AS character(255)),
              CAST(license_type AS character(255)),
              CAST(url AS character(255))
          FROM SDMI_Inventory" \
    sdmi_inv.gdb \
    WFS:"http://wfs.alaskamapped.org/inventory"

This worked great until 10.0 Service Pack 4, at which point ArcMap / ArcCatalog started reporting errors similar to this when you try to add / preview the layer:

Error opening feature class
The item does not have a definition. [SDMI_Inventory]
The item does not have a definition.
The XML being loaded could not be parsed. [physical name: SDMI_INVENTORY]
The item does not have a definition.
The XML being loaded could not be parsed. [physical name: SDMI_INVENTORY]

I've got the latest (1.2) API and the latest GDAL (1.9.0), x86_64 GNU/Linux.

Anyone have any ideas how to get this to work again?


texts: Calculating closest distance from location to given route?

texts:
I want to calculate the closest distance from a location (lat/lon) to a given route. The route consists of approximately 20.000 location points.

I could calculate the distance to every location point and check the results for shortest distance, but that seems quite CPU-intensive as I receive a new location every second.

Does anybody know an algorithm that is covering my problem?


texts:
We have to develop an application which will have lot of rich ui along with a map. For the map, I am using the arcgis javascript api which includes dojo. For UI components other than the map, I want to use a copy of dojo which I will host locally. How can I load both of these?
   I read this blog esri dojo, but this works with api 1.6, and I am using 3.6.


texts: What causes this unexpectedly poor GeoServer/PostGIS/OpenLayers performance?

texts:
I am having a tough time making this map run smoothly (play around with it for a bit, you'll see). I have taken a few tips from OSGeo's Making Maps Fast document, but I think there's another issue because:

A) sometimes the map runs just fine, but when I refresh the page it takes a very long time to load and/or times out on some layers, and

B) occasionally when zooming in/out, the Bing basemaps change appropriately but my layers stay at their original scale.

I have not yet taken the step of simplifying/generalizing geometries because, although it may help, I don't think that's the issue. I'm only bringing in a few layers on the initial load, and they're not overly complex, at least not complex enough to slow things down this much IMO. Also, I briefly employed GeoWebCache, but that didn't make much difference either.

In terms of additional performance improvement techniques, I have also followed GeoServer's container considerations and my JAVA_OPTS look like so:

grep -i useparallelgc /etc/*/*
/etc/init.d/tomcat6:export JAVA_OPTS="-server -Xmx256M -Xms48m -XX:SoftRefLRUPolicyMSPerMB=36000 -XX:MaxPermSize=128m -XX:+UseParallelGC"
/etc/java/java.conf:export JAVA_OPTS="-server -Xmx256M -Xms48m -XX:SoftRefLRUPolicyMSPerMB=36000 -XX:MaxPermSize=128m -XX:+UseParallelGC"

Here's my environment:

GeoServer 2.1.1
PostGIS 1.5.3
PostgreSQL 8.4ish
OpenLayers 2.11
GeoExt 1.1
ExtJS 3.4.0
Amazon EC2 Linux Micro Instance (awww, 613 MB of memory)

If it helps, here is my GeoServer server status:

Data directory      /usr/share/tomcat6/webapps/geoserver/data
Locks           0
Connections     2
Memory Usage        144 MB
JVM Version     Sun Microsystems Inc.: 1.6.0_29 (Java HotSpot(TM) 64-Bit Server VM)
Available Fonts     GeoServer can access 97 different fonts. Full list of available fonts
Native JAI      true
Native JAI ImageIO  false
JAI Maximum Memory  113 MB
JAI Memory Usage    0 KB
JAI Memory Threshold    75.0
# of JAI Tile Threads   7

JAI Tile Thread Priority        5
ThreadPoolExecutor Core Pool Size   5
ThreadPoolExecutor Max Pool Size    10
ThreadPoolExecutor Keep Alive Time (ms) 30000
Update Sequence             236

You can look at my code on the website itself, but here is an example of the parameters I'm using for all my WMS layers (tried them as tiled layers, which didn't seem to make much difference):

var wildSingleTile = new OpenLayers.Layer.WMS(
    "Wilderness Areas",
    "http://50.112.123.161:8080/geoserver/OR_wild_wkspc/wms",
    {
        LAYERS: 'OR_wild_wkspc:OR_Wilderness',
        STYLES: '',
        transparent: true,
        format: 'image/png8'
    },{singleTile: true, isBaseLayer: false}
);

I realize that Amazon's micro instance is not the speediest thing in the world, but I don't think it should be THIS slow.

I am quite the noob at just about everything here, including JavaScript and Linux, so any tips and suggestions would be appreciated!

P.S. Sorry for the hideous layer colors, it's not final...


texts: Checking which layer edited under specific version?

texts:
I have approx 100 feature class in SDE database. Out of this someone edited 10 in "test_1" version feature class(assumed). Now i want to know with script or manual which feature class is edited under "test_1" version.

How can I check this?


texts:
At present I have a dataset which is broken down by Electoral Divisions (circa 3,500). However, due to issues with some of my variables I need to choose a larger spatial unit - Local Electoral Area (LEA). An LEA is an aggregation of a number of different Electoral Divisions. The issue with the LEA shapefile, is the attribute table does not provide a lot of information of the region.

What I would like to do is aggregate my various Electoral Division variables to an LEA level. However, I believe this will be extremely difficult. Instead what I would like to do is join by ED's by LEA so I'll know what ED is in what LEA. The issue here is at present I have no identifier which is common between both shapefiles.

I originally was going to geocode by Electoral Divisions and than merge them to the LEA shapefile, however, 3,500 is quiet an amount, as well as ED names are quiet common with other areas in the world.

Is there some way to overlay the two shapefiles so my ED's are imbedded within the LEAs. Both shapefiles are working of the same coordinate system, so I imagine this is possible?

I'm using QGIS


texts: How can I keep my GEOID a 12 character string when adding csv layer in QGIS?

texts: I am trying to join my csv data to TIGER census block group shapefiles. Every time I bring the table into QGIS it changes the string for the GEOID into a number so the join comes up null. My csv only has the block group ID to make the link, no address info or coordinates.


texts: Maximum Liklihood Image Classification; Does this algorithm consider RGB Band values only?

texts: if it is correct that Maximum Liklihood classification algorithm takes into account RGB values of the bands only, then is it useless to stack other bandslike NIR SWIR etc prior to go ahead and use stacked image to retrieve landcover classes?


texts:
I created on a Wordpress site for my final thesis and need to add a map. I do not even know where to start.  I am so confused.  I do not need to do anything fancy (yet)  I just want to map about 30 spots in the South Bronx.
I was using leaflet but I was told that Carto was the way to go.

I found this https://jsmit8.carto.com/tables/breweries/public/map
It has current breweries.  I want to create a map with the breweries that no longer exist.


texts: ArcGis 10 interpolate Z between polyline's vertices (IZ.InterpolateZsBetween ?)

texts:
Here's what's happening:

I want to create a 3D line using only the height of some GSP points. I need to be able to extrapolate the vertice's height between where it's known to the one where it's unkonwn.

What I have tried:

Creating a TIN from GPS points: where the line meanders and there is no GPS points, the height is not correct.
Linear Referencing tools: but as originally the Polyline Z has no value interpolated between each vertex, it doesn't seem to do what is needed
Conversion of Point to raster at a 1 m cell (Nodata outside the cell where the GPS points lay) and extraction of the Z value to the vertex where they overlap the cell with value, leaving the not-overlapping vertices without Z value.

Is this what I need?:

I think I should use something like "IZ.InterpolateZsBetween" but I am useless at coding.
Could someone please explain me how to implement this fonction?

Many thanks


texts: Is it possible to install the PostGIS to the PostgreSQL (binary version)?

texts:
I'm using the PostgreSQL 9.2.4 which is the no installation version.I run this by using the command "Program Files/PostgreSQL/9.0/bin/pg_ctl" -D "Program Files/PostgreSQL/9.0/data" -l "Program Files/PostgreSQL/9.0/log/pgsql.log" start"!

And I tried to install the PostGIS 1.5.5 into the PostgreSQL server. It is not work! The template_postgis is created, but no tables was generated in the table folder!

Would like to confirm, does the PostGIS need to be installed to installion version instead of binary version? Please advise.

Thanks.


texts: hope anyone could give me some hints here. I tried to download free MODIS data from their page(resolution 500m, 1km). Im interested in NDVI and temperature data. But when i load the hdf files into my project file I just get weird layers :/.   This for example should be the temperature layer, but that's nothing I can work with. You have a clue if these are the right files, or if not, if there are better ones?


texts: ArcScene Export to 3D Webscene polygon issue?

texts:
I have been putting together a 3D basemap for our campus, Washington State University, and have been trying to share my model using the export to 3d web scene arcToolbox tool.  For the most part is has worked pretty well except for my StreetPavement feature class (polygon).  For some reason, blank spaces in between polygons are getting filled.  In other words, the areas that would be filled by my LandscapeArea polygon class are being filled with gray polygons.  I have attached a screenshot to help illustrate.  If I turn off the StreetPavement layer in the web scene viewer it looks fine.  The strange this is that the LandscapeArea class doesn't experience this same phenomenon.  I have tried about twenty different things to the StreetPavement class, but always end up with the same result.  For example:

Merging all of the LandscapeArea class, then performing a clip on the StreetPavement
Creating a new featureclass, then coping the features from the StreetPavement into it.
Creating a raster layer from the StreetPavement class then displaying it in the model.  (This sort of works but doesn't look that great)

I have even gone so far as to copy the features within the StreetPavement class and copying them into the LandscapeArea class.  The result ends up looking the same.  So I know the issue has to do something with the polygons in the StreetPavement class, but I have no idea what it is or what else to try.

Dark gray area filling in between polygons in the StreetPavement Class.

Without the StreetPavement Layer turned on.

Everything looks great in ArcScene.


texts: Optimizing Arcpy Code to Generate Points Along a Line

texts:
I have a line layer that shows roads across the country and I also have a table with road ID numbers that represent where potential stoplights would be placed. I have been able to use the shared ID numbers between the two feature classes to generate a point at the middle of a road segment and then add that point to a larger point layer. I am then able to calculate the geometry of the point layer and get lat long coordinates, which are essential to the analysis that I am doing. Here is what I have so far:

import arcpy, math, os, csv
arcpy.env.overwriteOutput = True
arcpy.env.workspace = "C:\Users\dsala\Desktop\GISDATA\SE_Br.gdb"

outTBL = 'ComID_Info'
in_Tbl = 'C:\Users\dsala\Desktop\GISDATA\SE_Br.gdb\light_ComID'

Selection_field_polyline = "COMID"
Selection_field_table = "V2_ComID"
list_Lights_COMID = []
list_Road_COMID = []
cur_Tbl = arcpy.SearchCursor(in_Tbl)

for num in cur_Tbl:
    list_Plant_COMID.append(num.getValue(Selection_field_table))
    polyline_layer = arcpy.MakeFeatureLayer_management('Road_network', "temp_Road_network")

cur_line = arcpy.SearchCursor(polyline_layer)
for row in cur_line:
    list_River_COMID.append(row.getValue(Selection_field_polyline))

S1 = set(list_Plant_COMID)
S2 = set(list_River_COMID)
intersect = S1.intersection(S2)
arcpy.CreateFeatureclass_management('C:\Users\dsala\Desktop\GISDATA','Lights_pt', 'POINT', 'temp_points_layer')

for x in intersect:
layer1 = arcpy.management.SelectLayerByAttribute( polyline_layer, "NEW_SELECTION", "\"COMID\" = {}".format(x))
Selection = arcpy.MakeFeatureLayer_management(layer1, "temp_road_location")
arcpy.GeneratePointsAlongLines_management(Selection, 'temp_points_layer', 'PERCENTAGE',  Percentage=50)

temp_pts_scur = arcpy.SearchCursor('temp_points_layer')
facilities_pts_icur = arcpy.InsertCursor('Lights_pts')
for pt in temp_pts_scur:
    facilities_pts_icur.insertRow(pt)

arcpy.AddGeometryAttributes_management ('Lights_pts', 'POINT_X_Y_Z_M')

I have been testing this code out on a few states, however when I expanded it nationwide it would never finish. It would spend too long trying to redraw the national road line layer every time my code selected a new line segment and would eventually crash. I am looking for a way to make this process more efficient. Ideally, there would be a workaround where I would not need to use the 'SelectLayerbyAttribute' Tool. I can't think of way to generate a new point layer that wouldn't involve selecting features to then use the 'GeneratePointsAlongaLine' tool.


texts: leaflet / omnivore - display omnivore CSV markers by property value else no display

texts:
I asked a question earlier about displaying different points with different symbols based on their values for a specific property. With the help of another user I was able to create if-else statements to the effect of--

window.dataWR = omnivore.csv('myPoints.csv')
    .on('ready', function(layer) {
    markerArray = [];
    this.eachLayer(function(marker) {
        if (marker.toGeoJSON().properties.Category == "Type A") {
            marker.setIcon(L.icon({
                iconSize: [18, 18],
                iconAnchor: [9, 27],
                popupAnchor: [1, -24],
                iconUrl: 'c2.png'
            }));

        }
        else if (marker.toGeoJSON().properties.Category == "Type B" ) {
            marker.setIcon(L.icon({
                iconSize: [18, 18],
                iconAnchor: [9, 27],
                popupAnchor: [1, -24],
                iconUrl: 'c3.png'
            }));

        }
        else {
            marker.setIcon(L.icon({
                iconSize: [0],
                iconAnchor: [0],
                popupAnchor: [0],
                iconUrl: 'c1.png'
            }))}

      return marker.bindPopup(marker.toGeoJSON().properties.Category + ', ' +
                marker.toGeoJSON().properties.Location);
    });
});

The one problem I have with this setup however is that the else conditional kind of "cheats" in that it does display even though the symbol size is set to zero. So if I have multiple statements like this, I'm always rendering the whole dataset even if the user can't see it, slowing down my map considerably. Is there a way of accomplishing this so only those data with that categorical value are displayed rather than the others having 0px symbols? Something to the effect of else "return false" or else "display:none"?


texts:
I want to use ogr2ogr command line to convert KML files to shapefile file, however, I can't get correct shapefile when I type the ogr2ogr convert command.

However, I can succeed in converting the continent.shp to correct KML file inside ArcGIS software by means of ogr2ogr. Besides, the conversion between acquired KML file in last step to shapefile.

I guess maybe my KML file is not the strict standard format. So I paste the KML file contents.

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
  xmlns:gml="http://www.opengis.net/gml"
  xmlns:xfdu="urn:ccsds:schema:xfdu:1"
  xmlns:safe="http://www.esa.int/safe/sentinel-1.0"
  xmlns:s1="http://www.esa.int/safe/sentinel-1.0/sentinel-1"
  xmlns:s1sar="http://www.esa.int/safe/sentinel-1.0/sentinel-1/sar"
  xmlns:s1sarl1="http://www.esa.int/safe/sentinel-1.0/sentinel-1/sar/level-1"
  xmlns:s1sarl2="http://www.esa.int/safe/sentinel-1.0/sentinel-1/sar/level-2"
  xmlns:gx="http://www.google.com/kml/ext/2.2">
  <Document>
    <name>Sentinel-1 Map Overlay</name>
    <Folder>
      <name>Sentinel-1 Scene Overlay</name>
        <GroundOverlay>
          <name>Sentinel-1 Image Overlay</name>
          <Icon>
            <href>quick-look.png</href>
          </Icon>
          <gx:LatLonQuad>
            <coordinates>-113.892487,68.977577 -124.235062,69.902451 -122.225281,73.452110 -109.974060,72.393707</coordinates>
          </gx:LatLonQuad>
       </GroundOverlay>
    </Folder>
  </Document>
</kml>

I don't know the reason to the conversion problems. Could you give me some advice?


texts: How to resample raster exportation with ArcObjects

texts:
I noticed the resolution of a programmatically exported raster was not as good as the same image exported with the ArcMap GUI. I set the resample ratio in the GUI to 3 and the quality seemed identical. I tried to do the same thing in code but it doesn't seem to be changing the quality at all:

private static IExport GetExporterObject()
{
    IExport newExporter = new ExportPDFClass();
    IOutputRasterSettings rasterSettings = (newExporter as IOutputRasterSettings);
    rasterSettings.ResampleRatio = 1;
    return (IExport)rasterSettings;
}

private static IEnvelope GetExporterRectangle(double ratio, IActiveView activeView, out tagRECT exporterRectangle)
{
    IEnvelope pixelBox = new EnvelopeClass();
    pixelBox.XMin = activeView.ExportFrame.left * ratio;
    pixelBox.XMax = activeView.ExportFrame.right * ratio;
    pixelBox.YMin = activeView.ExportFrame.top * ratio;
    pixelBox.YMax = activeView.ExportFrame.bottom * ratio;

    //I did the following code this way because I found that converting the ratio to int
    //was rounding it down before multiplying and giving me the wrong dimensions
    exporterRectangle.left = (int)(activeView.ExportFrame.left * ratio);
    exporterRectangle.right = (int)(activeView.ExportFrame.right * ratio);
    exporterRectangle.top = (int)(activeView.ExportFrame.top * ratio);
    exporterRectangle.bottom = (int)(activeView.ExportFrame.bottom * ratio);

    return pixelBox;
}

public static void ExportMapPage(IMxDocument mxd, MapPage mapPage, string outputPath)
{
    mxd.ActiveView = mxd.PageLayout as IActiveView;
    IActiveView activeView = mxd.ActiveView;
    activeView.FocusMap.MapScale = mapPage.Scale;
    mxd.ActiveView.PartialRefresh(esriViewDrawPhase.esriViewAll, null, activeView.Extent);
    string output = String.Format("{0}\\{1}{2}", outputPath, mapPage.Year, GetFileExtension(mapPage.FileFormat));

    IExport exporter = GetExporterObject();
    exporter.ExportFileName = output;
    exporter.Resolution = mapPage.DPI;

    tagRECT exporterRectangle;
    exporter.PixelBounds = GetExporterRectangle(mapPage.Ratio, activeView, out exporterRectangle);

    int hdc = exporter.StartExporting();
    activeView.Output(hdc, (int)exporter.Resolution, ref exporterRectangle, null, null);
    exporter.FinishExporting();
    exporter.Cleanup();
    mxd.UpdateContents();
}

I also noticed in the documentation there appears to be a discrepency.

This link says DisplayTransformation is the only coclass that implements IOutputRasterSettings, while this link says many do. Does anyone know which is correct? I'm using ArcGIS 10.1 SP1 and SDK 10.1 but I'm not sure how that reflects the ArcObjects version I'm using.

EDIT: It appears that DisplayTransformation is the coclass to use. I found that casting IActiveView.ScreenDisplay.DisplayTransformation to IOutputRasterSettings told me that the ResampleRatio property was in fact set to 3 by default. However even when I added this code immediately before the exporter.StartExporting() line, the output still looks the same as the ratio being set to 3.

IDisplayTransformation dispTrans = activeView.ScreenDisplay.DisplayTransformation;
(dispTrans as IOutputRasterSettings).ResampleRatio = 1;

In the documentation it specifically says at the bottom that changing this property will affect the results of the IActiveView.Output method, so I assume there is something wrong with my code.


texts: How can I register a table from an external source with a geodatabase so that I can create relationship classes?

texts:
I have two databases. The first is an enterprise geodatabase. The second is a SQL Server database that I have read only permissions to on another server. I am using ArcMap and ArcServer 10.4.1.

I am able to establish database connections to both database servers and create relates and joins without issue.

My goal is to use relationship classes with data I publish to services that will end up in ArcGIS portal. I know I have to use relationship classes since joins and relates won't work with portal.

I was under the impression that I could register tables from my external database with my geodatabase with the "Registerwithgeodatabase" tool. Then, I planned on creating my relationship classes.

When I run the tool I get "001400 : Only the data owner may execute register with geodatabase". Since there is only one input, it seems as though the tool may be trying to register the table from my external database with a non-existent geodatabase on the same server.

What am I doing wrong here and how can I register an external database table where I have readonly permissions with my geodatabase so that I can create relationship classes?

Edit: My external database contains only non-spacial data.


texts: Using two map frames on one layout in ArcGIS Pro with ArcPy

texts:
I want to create a map series that exports two different sheets on one page. For example Page 1 would have sheet 1 & 2 and Page 2 would have sheet 3 & 4. I used to do this with python in ArcMap, but it seems the method to do this has changed. The way it is now only uses 1 frame.

Below is what I managed so far:

import time
import arcpy, os

aprx = arcpy.mp.ArcGISProject("CURRENT")
lyt = aprx.listLayouts("13 Dual Panel Route Map")[0]
m = aprx.listMaps("ROUTE MAP")[0]
ms = lyt.mapSeries
df1 = lyt.listElements("MAPFRAME_ELEMENT","Map Frame 01")[0]  #MAKE SURE LAYER NAME MATCHES DATAFRAME
df2 = lyt.listElements("MAPFRAME_ELEMENT","Map Frame 02")[0] #MAKE SURE LAYER NAME MATCHES DATAFRAME
lyr = m.listLayers("MapIndex")[0]
index = lyr # my DDP index layer

# ColumnNumber          0      1      2
features = {row[0]:[row[1],row[2],row[3]] for row in arcpy.da.SearchCursor
     (index,['Id','SHAPE@','Scale','Name'])} # all index features(EDIT AS NECESSARY)
outpath = "A:\PROJECTS\WABASH VALLEY POWER\Enterprise South - Anson North\\05 PDF\\" #********CHANGE LINE HERE****************

#This part runs through the DDP, changing the pages and exporting them.
for pageNum in range(1,ms.pageCount,2):
    ms.currentPageNumber = pageNum
    zoom = "Id = 1+1"
    arcpy.SelectLayerByAttribute_management(lyr,"NEW_SELECTION",zoom)
    df2.camera.setExtent(df2.getLayerExtent(lyr, True, True))
    ex_output = outpath + features[pageNum][2] + ".pdf" #EDIT AS NECESSARY
    #df1.scale = features[pageNum][1]
    #df2.extent = features[pageNum+1][0].extent
    #df2.scale = features[pageNum+1][1]
    print ("Exporting " + features[pageNum][2])
 #(out_pdf, {resolution}, {image_quality}, {compress_vector_graphics}, {image_compression}, {embed_fonts}, {layers_attributes}, {georef_info}, {jpeg_compression_quality}, {clip_to_elements}, {output_as_image}#
    lyt.exportToPDF(ex_output,300,'BEST',False,'NONE',True,'LAYERS_ONLY',True,False,False,False)
    time.sleep(0)
del aprx


texts: Can Leaflet ZXY tiles be used in QGIS 3?

texts:
I have created tiles from a drone flight, loaded them onto Amazon S3 storage and used them in a leaflet web page served from an AWS EC2 instance.    That is all working fine.  I would my desktop QGIS 3.10 application like to use the same tiles stored on the S3 bucket in ZXY folders.

I tried to ust the XYZ Tiles point to the directory URL similar to this => https://my-maps.s3.us-east-2.amazonaws.com/MAP1/2020/{z}/{x}/{y}.png.   This does not work because I probably need a tile server running there.

I do not really want to have 2 copies of the map on different sites if I can do it from one spot.  Is there a way for QGIS 3 to use these tiles from Amazon S3 bucket storage?

Reference 1: https://gis.stackexchange.com/questions/304483/is-it-possible-to-serve-vector-tiles-directly-from-s3

Reference 2: https://github.com/addresscloud/serverless-tiles


texts: Fading symbology at edges using ArcGIS Desktop?

texts:
A general problem for me with digitized material is that it can be too precise. Sometimes I just need to show a line or surface without sharp borders. I don't regularly need to visualize the exact extent of, say, a wetland to meter or dm level, only show that it exists in approximately that area.

Now I have lines (that I also have buffered 50 m) which I need to "smooth" at the edges. I want the center of the line/buffer to be strong in color. Further away from the center it should be more transparent, and eventually fade away.

How can I accomplish this?


texts:
Does anyone have experience using QGIS to submit shapefiles to the Minnesota DOT (MNDOT)?

I am using EPSG:26915 but got rejected due to an error in the CRS.  MNDOT Specs call for "UTM coordinates, Zone 15E extended, NAD83, GRS80, horizontal units = meters, vertical units (if used) = feet."


texts: Getting full path of layer selected in drop-down box for ArcGIS Python script tool?

texts:
I have a Python script-based tool in ArcGIS. One of the parameters for the tool is an input raster file, which I've set (using the properties of the tool within the toolbox) to be a Raster Layer. This means that when I execute the tool I get the option to select that parameter from a dropdown list of currently loaded raster layers, or by navigating to find a new raster layer within the filesystem.

I am then using the standard arcpy.GetParameterAsText(0) code to get the parameter into my Python script. This works fine when I select a file by browsing through the filesystem, as the text in the dialog is the full path to the file, but when I just select from the dropdown list the text I get back is just the filename (eg. file.tif).

The code I am running needs to know the full path to the file - how do I get this?

If I could assume that the file was always in the workspace then I could append the filename to arcpy.env.workspace, but I can't assume that. Do I need to iterate through all of the layers that are loaded until I find one with the same name, and then find its full path, or is there an easier way?


texts:
I need only basic osm2pgsql usage. Using postgres user and sandbox  (created) database.

   sudo -u postgres osm2pgsql -s -U postgres -d sandbox brazil-latest.osm.pbf

NOTES

Using UBUNTU 14 LTS, osm2pgsql --version say "osm2pgsql SVN version 0.82.0 (64bit id space)". At SQL terminal (psql) the select version() says "9.3.9 on x86_64-unknown-linux-gnu".

Prepare, all steps:

wget -c http://download.geofabrik.de/south-america/brazil-latest.osm.pbf

wget -c http://download.geofabrik.de/south-america/brazil-latest.osm.pbf.md5

md5sum -c brazil-latest.osm.pbf.md5
# OK

sudo -u postgres psql sandbox
# ...     SELECT PostGIS_version();  -- 2.1 OK
\q

sudo -u postgres osm2pgsql -s -U postgres -d sandbox brazil-latest.osm.pbf

last command errors:

osm2pgsql SVN version 0.82.0 (64bit id space)

Using projection SRS 900913 (Spherical Mercator)
Setting up table: planet_osm_point
NOTICE:  table "planet_osm_point" does not exist, skipping
NOTICE:  table "planet_osm_point_tmp" does not exist, skipping
Setting up table: planet_osm_line
NOTICE:  table "planet_osm_line" does not exist, skipping
NOTICE:  table "planet_osm_line_tmp" does not exist, skipping
Setting up table: planet_osm_polygon
NOTICE:  table "planet_osm_polygon" does not exist, skipping
NOTICE:  table "planet_osm_polygon_tmp" does not exist, skipping
Setting up table: planet_osm_roads
NOTICE:  table "planet_osm_roads" does not exist, skipping
NOTICE:  table "planet_osm_roads_tmp" does not exist, skipping
Allocating memory for dense node cache
Out of memory for node cache dense index, try using "--cache-strategy sparse" instead
Error occurred, cleaning up

Try also sudo -u postgres osm2pgsql --create --cache-strategy sparse -s -U postgres -d sandbox brazil-latest.osm.pbf  but some error: "Allocating memory for sparse node cache
Out of memory for sparse node cache, reduce --cache size
Error occurred, cleaning up
"

Try also osm2pgsql -c -d sandbox -U postgres -H localhost -S brazil-latest.osm.pbf  that say "Usage error".

Try also osm2pgsql -U postgres --slim -C 20000 -d sandbox --host localhost --number-processes 20 ./brazil-latest.osm.pbf
but "Error: Connection to database failed: fe_sendauth: no password supplied".

sudo -u postgres osm2pgsql -U postgres --slim -C 20000 -d sandbox --host localhost --number-processes 20 ./brazil-latest.osm.pbf  say "Error: Connection to database failed: fe_sendauth: no password supplied".


texts: Apply reduction without time pixels reducing

texts:
I'm having a little bit of trouble reducing the NDVI band of my data, without reducing the time bands. What I want is to apply a 95th percentile reducer to the NDVI band of my image from Colorado and have the time stamp from the pixels that were selected for the composite image to carry over, but instead I am reducing BOTH the NDVI and system:time_start band.

//ADD QUALITY BANDS TO MAKE GREENEST PIXEL MOSAIC FOR LANDSAT 5
//This function masks clouds in Landsat 5 imagery.
var maskClouds = function(image) {
  var scored = ee.Algorithms.Landsat.simpleCloudScore(image);
  return image.updateMask(scored.select(['cloud']).lt(20));
};
//This function masks clouds and adds quality bands to Landsat 8 images.
var addQualityBand = function(image) {
  return maskClouds(image)
    //NDVI
    .addBands(image.normalizedDifference(['B4', 'B3']))
    //time in days
    .addBands(image.metadata('system:time_start'));
};

var addTime = function(image) {
  return image.addBands(image.metadata('system:time_start'));
};

// Map the function over the collection and display the result.
////print(collection.map(addTime));

//Map.addLayer(fldas_colorado_1990, {}, "max temp and precip data for 1990");
var collection_1990= ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA')
  .filterDate('1990-1-1','1990-12-31')
  .map(addQualityBand);
var only_NDVI_1990 = collection_1990.select('nd','system:time_start');
var ninety_five_percent_1990 = only_NDVI_1990.reduce(ee.Reducer.percentile([95])).clip(colorado);
Map.addLayer(ninety_five_percent_1990, imageVisParam, '1990 95');

Here is the link to my script.
https://code.earthengine.google.com/2664f992629c656230c8a9b531644660

Very basically, I want to be able to reduce bands in images, without reducing the system:time_start band, but still have the data from the reduced pixels carry over. Is this possible?

Thanks so much.


texts: how activate qgis features selection toolbar in the main window

texts: How can i activate and desactivate qgis features selection toolbar in the main window using a signal in a python plugin?


texts: Using python parser of ArcGIS field calculator to calculate field using information from subsequent row?

texts:
When using the field calculator, how do I get data from the subsequent row in  python?

Let's say I want to add a new field called head_diff. I want to take the current row's value for head, minus the subsequent row's value for head to get a head difference. Example for previous row:

for index, row in enumerate(rows,start=1):
    head_diff= !head! - row[index-1]

Someone suggested da.UpdateCursor, which looks promising. Will comment with procedure if I find a solution.


texts: Find all the values of a specific field of a layer based on location then select just those values to see where else they appear in the layer

texts: I have a layer of property ownership in a county and need to see all the property owners in a parish [parish value already defined to feature in the layer] then see what other property those owners in the parish own in the county. I can use QGIS 2.12 or Arc 10.3.1


texts: Creating duplicate features based on many-to-one conversion of related table?

texts:
I need to create duplicate features (in this case parcels), using a related table to populate ID fields.  In essence, I need to use a related table of parcel owners to create duplicate parcels, each with its own unique ID from the related owner table.  The related table already contains the many-to-one link to parcels, I just want to force a one-to-one relationship between the owners and parcels, by creating a duplicate parcel for each owner record.

I'm using ArcGIS Desktop 9.3.1.


texts:
I have a simple Leaflet map and would like to add a GPX file coming from a Garmin device. I use omnivore plugin and by default it gets a line. How could I represent points in each coordinate instead of a line?
var map = L.map('map').setView([41.39742, 2.16328], 13);

L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
        maxZoom: 25,
        attribution:'Map data &copy; <ahref="http://openstreetmap.org">OpenStreetMap'
    }).addTo(map);

var gpxURL = './gpx/Track_PESCA1-19 071255.gpx';
var gpx = omnivore.gpx(gpxURL).addTo(map);


texts: Map is not visible at initialization using react-leaflet

texts:
I have this react component:

import React, { Fragment } from 'react';
import L from 'leaflet';
import { Map, TileLayer, Marker, Popup } from 'react-leaflet';
import 'leaflet/dist/leaflet.css';

const customMarker = new L.icon({
    iconUrl: '/images/poemEditorTools/location-pointer-new.svg',
    iconSize: [56, 72],
    iconAnchor: [26, 72],
});

const MyPopupMarker = ({ content, position }) => (
    <Marker position={position} icon={customMarker} >
      <Popup>{content}</Popup>
    </Marker>
)

const MyMarkersList = ({ markers }) => {
    const items = markers.map(({ key, ...props }) => (
        <MyPopupMarker key={key} {...props} />
    ))
    return <Fragment>{items}</Fragment>
}

const markers = [
    { key: 'marker1', position: [51.5, -0.1], content: 'My first popup' },
    { key: 'marker2', position: [51.51, -0.1], content: 'My second popup' },
    { key: 'marker3', position: [51.49, -0.05], content: 'My third popup' },
]

class MapWithPoems extends BaseComponent {
    render() {
        return (
            <Map center={[51.505, -0.09]} zoom={13} style={{ height: "500px", width: "100%" }} >
                <TileLayer
                    url={
                        "https://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer/tile/{z}/{y}/{x}"
                    }
                />
                <MyMarkersList markers={markers} />
            </Map>
        );
    }
}

export default MapWithPoems;

The problem is when the page is downloaded the map is shrinket to the left side of the container. I have to resize the window to make the map to show up in the full width. See images.

What is wrong here?


texts: getSelectedExtent not working on joined layers

texts: When joining layers the getSelectedExtent function in ArcPy fails to get the extent.  Unjoining the layers and running it yeilds the expected result.  The joined layers are all located on an SDE.  Is this a known issue or is there something I'm missing?


texts: Error Dictionary.get: Dictionary does not contain key: NDVI

texts:
I am trying export NDVI result to CSV file for 35 shapefiles farmers in Earth Engine. But it did not work. I tried find any solution from internet. How can I solve this problem?

var Akdarya = ee.FeatureCollection("users/xodjayevshovkat/Akdarya_GEE_35");

//Define a geometry.
Map.addLayer(Akdarya);

// Make the geometry a feature collection.

// Import S2 image collection.
var S2 = ee.ImageCollection('COPERNICUS/S2')
  // Filter collection.
  .filterDate('2018-05-01','2018-09-30')
  .filterBounds(Akdarya)
  .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 20);

// Define a cloud masking function.
var maskcloud1 = function(image) {
  var QA60 = image.select(['QA60']);
  return image.updateMask(QA60.lt(1))
};

//

// Define an index function (return only NDVI).
var NDVI = function(image) {
  return image.expression(
    '(NIR - RED) / (NIR + RED)',
    {
      'NIR': image.select('B8'),
      'RED': image.select('B4'),
    }).rename('NDVI').copyProperties(image, image.propertyNames());
};

// Apply the cloud mask function.
var NDVIcol = S2.map(maskcloud1)
 // Apply the NDVI transformation.
 .map(NDVI);

// Calculate the per-pixel median/mean NDVI for all images in collection.
var NDVImed = NDVIcol.mean();

// Define vis params.
var NDVIvis = {min: -0.001, max: 0.9, palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
    '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'};

//

// create a collection of the NDVI

var NDVIclip = NDVImed.clipToCollection(Akdarya)

// get the mean value for the region from each image
var ts = S2.map(function(image){
  var date = image.get("DATE_ACQUIRED")//.format("Y-M-D");
  var mean = image.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: Akdarya,
    scale: 30
  });

// and return a feature with 'null' geometry with properties (dictionary)
  return ee.Feature(null, {'mean': mean.get('NDVI'),
                            'date': date})
});

// Export a .csv table of date, mean NDVI for watershed
Export.table.toDrive({
  collection: ts,
  description: 'NDVI',
  fileFormat: 'CSV'
});

'''


texts:
I plan to use for positioning and tracking of areas worked for underwater mining.

I am trying to draw rings using Rectangles Ovals Digitizing plugin.

But there is very little info available.


texts: Joining fields in Processing Modeler of QGIS

texts:
I am trying to join field district layer with output A .

How can I do this?

I added Vector field input but i cannot see the output A field to make the join by fields .

I have attached an image to illustrate this case.


texts: Is there a way for apply point.distance on an array instead of doing double for loop?

texts:
I am using the code below to find a distance between one point and a line, but I have a grid of points, and I want to optimize it so I don't have to do double for loop. Is it possible to apply point.distance() on an array?
My goal is to get a 2D array of the distances. But it will be a 5000*5000 array.
This is a fully working example.

import shapely.geometry as geom
import numpy as np
nr = nz = 1000
r_arr=np.zeros((nr,nz))
z_arr=np.zeros((nr,nz))
temp = np.empty_like(z_arr)
start_r = 0.0e-9
end_r = 100.0e-9
start_z = 0.0e-9
end_z = 100.0e-9
linspr = np.linspace(start_r,end_r,nr)
linspz = np.linspace(start_z,end_z,nz)

for i in range(nr):
    r_arr[i]=linspr
for j in range(nz):
    z_arr[:,j]=linspz
coords = np.zeros((nm,2))
# you can use any random line for coord, my values come from another function but for simplicity I will use a straight line

coords[:,0] = end_r*np.ones(nr)/2.0
coords[:,1] = np.linspace(start_z,end_z,nz)

line = geom.LineString(coords)

#r_arr, z_arr are 2D array which have r,z for the points.
for i in range(nz):
    for j in range(nr):
        temp[i,j] = geom.Point(r_arr[i,j],z_arr[i,j]).distance(line)


texts: Maintaining sorted unique values in value table of ArcGIS Pro Python script tool validation?

texts:
I have created a Python script tool in ArcGIS Pro 1.4.1 which has tool validation.

The fourth parameter on my tool dialog, which is of type string with multiple values allowed, after choosing some values looks like this:

I would like this list of value to be sorted and contain only unique values.

I think the updateParameters function should be able to achieve that.  My updateParameters function appears below.

The code has an if statement to cater for the function being called before any values have been added.  Then what I think it should be doing is:

Read the values from the ValueTable object of parameter 3 into a string of values separated by semi-colons
Create a list of those values by splitting on semi-colons
Make a unique sorted list
Join that list back into a string and load that string into the value table object.

At this point I would expect my five unordered values to have become four values sorted alphabetically.

My code gives no error, but when I change the value of a parameter, triggering the updateParameters function, it leaves the list with the original five members still unordered.

def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed. This method is called whenever a parameter
    has been changed."""
    if self.params[3].value is not None:
        valueString = self.params[3].value.exportToString()
        # Sort and make values unique
        valueList = valueString.split(";")
        sortedList = sorted(set(valueList))
        self.params[3].value.loadFromString(";".join(sortedList))

Can anyone see what I am missing or may be misunderstanding?

I have assumed that the string is being exported from and read by the value table using semi-colons but do not know where I can check that.

The answer by @AlexTereshenkov has all but solved this for me.  The lingering issue which may or may not represent a new problem (that may need a new question) is an error that it seemed to bring in:

The value table appears to be updating perfectly but it is throwing ERROR 000800 and saying that one or both of the values in the picture are not in the choice list but they both are.

def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed. This method is called whenever a parameter
    has been changed."""
    if self.params[3].altered:
        self.params[3].value = sorted(list(set(self.params[3].value.exportToString().split(";"))))

I suspect that I'm seeing the error when Alex is not because I have spaces in my strings when he does not.  I'll experiment some more tomorrow.


texts:
This is scenario: I am tourist in a new place, say city, and I want to use  smart phone to get interesting places around me. To do that I search OSM (using overpass API) for nodes and ways with tag "tourism", and query returns bunch of places. If there a way to rank them? I saw questions about "popularity rank" in OSM, and I understand it's not OSM purpose. However as a tourist I do not necessary want to go to a popular place (to spend time in line), I want to go to place with maybe nice view, a story, place related to a famous person.

I was thinking I could rank places by number of tags (the more tags, the more info, more likely there is something interesting there), by number of edits, by number of nodes in OSM way (if somebody spent time entering all tags and defining way, maybe it's worth visiting), by presence of tag with URL. Any other ideas?


texts: Create a grid layer over existing raster layer using python modules

texts:
I have a  DEM raster layer and want to create a square grid of size 10x10 and align into my raster using python modules. I have tried to implement below code. But I'm not sure of what geometry to be passed for creating all_bbox parameter in the below code.

import shapely.wkt
from shapely.geometry import box

def generate_tiles():

  # width and height of the tiles to generate
  WIDTH, HEIGHT = 10, 10

  # from user arbitrary polygon, could be buffered
  **all_bbox = shapely.wkt.loads('POLYGON (())')**
  all_area = all_bbox.area

  coords = all_bbox.exterior.coords

  # reset x
  reset_lon = coords[0][0]
  # x
  lon = coords[0][0]
  # y
  lat = coords[0][1]

  covered = 0
  bboxes = []

  col_idx = 0
  row_idx = 0

  while covered < all_area:

      bbox = box(lon, lat, lon + WIDTH, lat + HEIGHT)

      intersection = all_bbox.intersection(bbox)

      covered += intersection.area
      bboxes.append(bbox)

      # catch the very last top row, the first column may already intersect
      if col_idx == 0 and intersection.area < bbox.area:
        pass

      # increase Y and reset X
      elif col_idx > 0 and intersection.area < bbox.area:

        # increase row
        lat += HEIGHT
        lon = reset_lon

        col_idx = 0
        row_idx += 1

        continue

      # just increment column index
      else:
        col_idx += 1

      lon += WIDTH

  return bboxes


texts:
I'm trying to understand how Arel work.

How can I do this PostGIS query with Arel.spatial ?

SELECT
    ST_asText(ST_ClosestPoint(track.path::geometry, pta.lonlat::geometry))  AS closet_point
    FROM tracks track, pois pta
    WHERE track.id = 2
      AND pta.id = 2

pta is a point like t.geography "lonlat", limit: {:srid=>4326, :type=>"st_point", :geographic=>true}
track.path is a linestring like t.geometry "path", limit: {:srid=>4326, :type=>"line_string"}


texts: ArcGIS Intersect Tool not splitting lines as expected?

texts:
I performed Intersect Analysis on one polygon shapefile and one polyline shapefile. I didn't specify any tolerances, and both SHPs were projected into the same UTM before performing the analysis.

Expected behavior: one line criss-crossing a polygon will always be split into multiple lines

Actual behavior: lines are only sometimes split

Evidence: # lines input file = 2663; # lines output file = 2958; so, some splitting is happening. Image below is only one FID highlighted, but clearly is 3 segments and should be 3 separate FIDs

Am I using the wrong tool?


texts: Getting Exent of Queried Results using ArcGIS Python API?

texts:
Im trying to use the ArcGIS Python API to query a few results and gather the extent, I can query the results fine but as soon as I try to extract the extent I retrieve empty results.

Using one of the samples servers, heres my code:

import arcgis.gis
gis = GIS()

t = 'https://sampleserver6.arcgisonline.com/arcgis/rest/services/DamageAssessment/FeatureServer/0'
layer = FeatureLayer(t, gis)

results = layer.query(where='1=1', return_extent_only=True)

After I print results I get the following

{"features": [], "fields": []}

Am I missing something? Seems all to easy up until this point.


texts: How do I calculate proximity for multiple points?

texts:
I have two layers, one with lots, one with subway stations (NYC). There are to scenarios:

Case 1. I want to rank the lots relative to the proximity to subways.

(1.a.) The closer to closest station, the better or
(1.b.) The closer to multiple stations, the better.

Case 2. I want to evaluate the area surrounding the lots, say 500m radius. I would calculate the distance from the center of each block to the nearest subway station(s). For each block, the distance will be multiplied by the number of residents. Then the sum of those values will be divided by the total population in the area.

How can I do this in QGIS or/and PostGis?


texts: Listing all feature datasets and classes from multiple geodatabase into CSV file?

texts:
I have a folder having multiple geodatabase and I am calling that folder “Multiplegdb”
There are 3 geodatabase in a folder which are as follows:

ABC.gdb
UVW.gdb
XYZ.gdb

Each geodatabase have several ‘Feature Dataset’ and each ‘Feature dataset’ have several ‘Feature Classes’.

My task:

I would like to list all the features datasets and feature classes in each geodatabase in a folder into CSV. I also would like to get the metadata information from each feature class into CSV file. The format of CSV file will be:
GDB Name, Feature Dataset Name, Feature Classes name, Summary, Descriptions, Credits, Date Modified, File path.


texts: ArcGIS JS API 3.x - Update the default selected feature highlight styling when a pop up displayed

texts:
I'm displaying a pop up when a feature is selected on feature layer. I'm also highlighting the polygon boundary with a red line but when the pop up appears the highlighter changes to aqua color which i'm assuming is the default behavior. The highlighter goes back to red as soon as the pop up is closed .

Without PopUp

Is there a way to override the default feature highlighter behavior when the Pop Up is displayed or just not show that graphic at all?


texts: How do I name a file with looping raster's name?

texts:
I'm using arcpy and i need to save a raster file with the looping raster's name inside the output's name, in other words, if the looping raster is named "20170823" it must save like "Example_20170823".
I need to find something similar to 'Example_%20170823%' that we have in arcgis.

import arcpy
import os
from arcpy import env
from arcpy.sa import *

rasters = arcpy.ListRasters("*", "TIF")
arcpy.env.workspace = 'Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL'

for tiff in rasters: #the loop

    arcpy.MakeRasterLayer_management(tiff, "Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\temp_rast1.tif", "", "", "3")
    arcpy.MakeRasterLayer_management(tiff, "Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\temp_rast2.tif", "", "", "4")

    arcpy.CopyRaster_management("Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\temp_rast1.tif", "copy_rast1", "","", "", "", "", "", "", "", "", "")
    arcpy.CopyRaster_management("Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\temp_rast2.tif", "copy_rast2", "","", "", "", "", "", "", "", "", "")

    Property1Result = arcpy.GetRasterProperties_management ("copy_rast1", "MINIMUM")
    Property2Result = arcpy.GetRasterProperties_management ("copy_rast2", "MINIMUM")

    Property1 = Property1Result.getOutput
    Property2 = Property2Result.getOutput

    Correcao1 = Minus ("copy_rast1", "Property1")
    Correcao2 = Minus ("copy_rast2", "Property2")

    Rastdiv1 = Minus ("Correcao2", "Correcao1")
    Rastdiv2 = Plus  ("Correcao2", "Correcao1")

    NDVI_tiff = Divide (Rastdiv1, Rastdiv2,'Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL_PL\\teste\\NDVI_%tiff%.tif') # here the name to replace %tiff%

The NDVI_%tiff% must take the information of the tiff's name to add in the final output file's name .


texts: Update image visualization parameters from slider UI elements in Earth Engine

texts:
I'm trying to change the visualization parameters (min and max) of NDVI images in an Earth Engine App with slider controls. The result is an animation of styled images.

The animation part works well, but in order to see the variance in the images, the user must be able to change the min and max according to their region of interest with a slider.

My attempt to update image styling is on line 46, but it does not work as intended. Here's the code:

Map.setCenter(-102.24954, 36.10286, 14);

var palette1 = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163',
'99B718', '74A901', '66A000', '529400', '3E8601', '207401',
'056201', '004C00','023B01', '012E01', '011D01', '011301']
var min1 = 0.2;
var max1 = 0.8;
var images = ee.ImageCollection("COPERNICUS/S2")
    .filterDate('2018-01-01','2030-12-31')
    .sort('system:time_start', false)
    .filterBounds(Map.getCenter());
var visNDVI = {min: min1, max: max1, palette: palette1};
var ndvi = images.map(function(image) {
  var result = image.normalizedDifference(['B8','B4']).rename("ndvi");
  return image.addBands(result);});
var imagesNDVI = ee.ImageCollection(ndvi).select('ndvi').limit(8);
var animation = require('users/gena/packages:animation');

var label = ui.Label('Set Min NDVI Value');
var slider = ui.Slider({
  min: 0, max: 1, step: 0.05, value: 0.3,
  style: {stretch: 'horizontal', width:'300px'},
  onChange: updateLayer
});

var label1 = ui.Label('Set Max NDVI Value');
var slider1 = ui.Slider({
  min: 0, max: 1, step: 0.05, value: 0.7,
  style: {stretch: 'horizontal', width:'300px'},
  onChange: updateLayer
});

var panel = ui.Panel({
      widgets: [label, slider, label1, slider1],
      layout: ui.Panel.Layout.flow('vertical'),
      style: {position: 'top-right',width: '295px'
      }
    });
Map.add(panel);

function updateLayer(value){
  var min1 = slider.getValue();
  var max1 = slider1.getValue();
  Map.layers().reset();
  var desired =     imagesNDVI.updateMask(imagesNDVI.get(visNDVI).get(min).set(min1).get(max).set(max1));
  Map.addLayer(desired, {palette: palette1});
}

animation.animate(imagesNDVI, {maxFrames: 50, vis: visNDVI, timeStep: 500});

Link to my GEE script


texts: adding non spatial attributes using geoserver

texts:
I am using Openlayers + geoext to perform transaction operations in postgis database using geoserver . i can insert spatial data into postgis but cannot insert other data like road name, road type in my road database.

the openlayers generates this code during transaction

<wfs:Transaction xmlns:wfs="http://www.opengis.net/wfs"
service="WFS"
version="1.1.0"
xsi:schemaLocation="http://www.opengis.net/wfs
http://schemas.opengis.net/wfs/1.1.0/wfs.xsd"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <wfs:Insert>
<feature:Complaint xmlns:feature="http://geoserver/postgis">
  <feature:geom>
    <gml:MultiCurve
    xmlns:gml="http://www.opengis.net/gml"
    srsName="EPSG:4326">
      <gml:curveMember>
        <gml:LineString>
          <gml:posList>85.292019896027 27.722567213659 85.303131631935 27.72224973549 85.303290371019 27.72478956084
          </gml:posList>
        </gml:LineString>
      </gml:curveMember>
    </gml:MultiCurve>
   </feature:geom>
  </feature:Complaint>
 </wfs:Insert>
</wfs:Transaction>

in geoserver demo i add fid and id so that it looks as

<wfs:Transaction xmlns:wfs="http://www.opengis.net/wfs"
service="WFS"
version="1.1.0"
xsi:schemaLocation="http://www.opengis.net/wfs
http://schemas.opengis.net/wfs/1.1.0/wfs.xsd"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <wfs:Insert>
<feature:Complaint xmlns:feature="http://geoserver/postgis">
  <feature:geom>
    <gml:MultiCurve
    xmlns:gml="http://www.opengis.net/gml"
    srsName="EPSG:4326">
      <gml:curveMember>
        <gml:LineString>
          <gml:posList>85.292019896027 27.722567213659 85.303131631935 27.72224973549 85.303290371019 27.72478956084
          </gml:posList>
        </gml:LineString>
      </gml:curveMember>
    </gml:MultiCurve>
   </feature:geom>
 **<id>650</id>
      <fid>400</fid>**
  </feature:Complaint>
 </wfs:Insert>
</wfs:Transaction>

and i get what i want. my question is how can i post the 2nd request (with id and fid) instead of first using openlayers and/or geoext?


texts:
I'm currently tasked with finding a way to deploy WFS services to a number of clients. The services could be consumed by both desktop GIS software e.g. MapInfo/ArcGIS or by Web GIS solutions.

Is it possible, for a WFS service, to define e.g. user name/password or token to pass along for verification of access rights? (top notch security is not necessary)

Have been looking in OGC WFS specification, but I can't seem to find any relevant info on this.


texts: How generate transparent tiles via gdal2tiles.py?

texts:
I have an ecw file and I want create tile via gdal2tiles.py. The file is as follow:

Note: Green color is style of a polygon layer that is under the ecw layer.
I want to transparent nodata sections in tile generating, How can I do it?


texts: Data from WFS: extract feature and not 3 band RGB (QGIS)

texts:
I am trying to extract the the data from the following WFS:

https://globalatlas2.masdar.ac.ae/geoserver/gwa/wfs

I am interested in the layer "Global Wind Atlas wind speed mean 1km at 100m height DTU 2015". I am using QGIS, version 2.18.15.

When I right click and choose "Save as" for the layer I can download it. However, the issue I am having is that it downloads a 3-band raster of pixel RGB values rather than the underlying feature. This is shown in the attached image using the "Identify Features" tool: the top image is the data that I want, from the WFS, but the bottom is actually what is saved.

How do I save the feature as it is defined in the WFS, and not the Pixel RGB values? I tried both the "Raw Data" and the "Rendered Image" save modes but both produced the same result.


texts: QGIS 3.2 how to save form attributes to a geopackage file?

texts: I have created a geopackage file in QGIS 3.2.  I have made amendments using the layer properties>form attributes to create alias names and change the widget types to checkboxes. When I send these files for remote workers to use it doesn't save any of these changes.  It only works if I send the project which I don't want to do.  Anyone any ideas on how to save these changes to the geopackage file.


texts:
Currently my script kicking out:

  Runtime Error: An expected Field was not found or could not be retrieved

I did some checks and the spelling of the fields are correct and the fields do exist in the feature layer.  The line below is the error.

for mat, sId in sCur:

listFields = [f.name for f in arcpy.ListFields("Target_Services")]
print listFields
for field in [srvMatFld, serviceID]:
    if field not in listFields:
        print field
with arcpy.da.SearchCursor ("Target_Services", [srvMatFld, serviceID]) as sCur:
    for mat, sId in sCur:


texts: Way to turn on a layer, and turn off the rest of the desired layers with a button or switch in Experience Builder?

texts:
I am attempting to publish a map in Experience Builder that has indoor floorplans of buildings. Would anyone know of a way to design a button or a switch that could toggle on a layer, and then turn off the other layers in question? I've seen some stuff for a floor switcher but that's in webapp builder.

Example: I have 1st floor, 2nd floor, and 3rd floor as layers. I'd like to design a widget similar to the one found on University of Oregon's indoor map, found here: https://map.uoregon.edu/[![enter image description here]2]2
Ideally it would look like that, but if not even a way to distinguish each floor , with buttons that look like 1 2 [3] for each floorplan. I'd like the user to have the ability to turn on a floor layer with a click, and that same click would then turn off the other floor layers.
Thanks in advance.


texts: Does ST_GeomFromText find exact matches for provided string or it finds similar, but not 100% same string?

texts:
SELECT * FROM centroids WHERE centroid = ST_GeomFromText('POINT (21.6198205 41.9418615)', 4326);
Record have centroid geometry equal to POINT (21.619820510769063 41.94186153740355)

Running on PostgreSQL with PostGIS 2.3 db shows me this record.
Running on PostgreSQL with PostGIS 2.4 db doesn't show me anything.

When I adjust query to exact geometry as in the record, so it is
SELECT * FROM centroids WHERE centroid = ST_GeomFromText('POINT (21.619820510769063 41.94186153740355)', 4326); it shows me record on both databases.
I couldn't find anything about how exactly does ST_GeomFromText work from technical side.
Is it correct behaviour?


texts: Getting the height and width of labels in QGIS?

texts:
Is there a way to get the height and width of labels in the field calculator? Perhaps with the help of the function editor? Units are in meters.

See also: Creating leader lines for moved labels


texts:
I am adding layers to my map dynamically based on whats available in the MapServer. I am looping over them so I can add them in the LayerControl.
My issue is that I dont have access to the layers outside of the scope of the function which I need to run a map.hasLayer() function. What can I do to access the layers outside the fetch function to perform removeLayer and other things?

let url = "https://sampleserver6.arcgisonline.com/arcgis/rest/services/Energy/Infrastructure/MapServer";
  let serviceUrl = url+'?f=pjson';
  fetch(serviceUrl, {method:"get"})
  .then(function(response) {
    return response.json();
  })
  .then(function(esriService) {
    esriService.layers.forEach(function(layer) {
      let layerName = layer.name.replace(/\s/g, "");
      let fullName = layer.name // return string
      layerName = L.esri.dynamicMapLayer({
        url: url,
        layers: [layer.id]
      }).addTo(map);
      layerControl.addOverlay(layerName, fullName)
    })
  });

map.on('click', function (e) {
  if (!map.hasLayer(mydefinedlayer)) {
    console.log("Layer not displayed on map")
  } else {
    console.log("Layer displayed on map")
  }
});


texts: How to connect to MSSQL-Server using Geoserver?

texts:
I have Geoserver on Ubuntu 13.04.  On another machine I have MS SQL-Server 2008 R2.  Now I want to connect it together. Google says that its possible to connect to SQL-Server from Ubuntu using ODBC but Geoserver works through JDBC.

In Windows its easy, you just to need to download MS JDBC driver. But what do you do in Ubuntu?


texts: merge raster images: error message procedure not found

texts:
I am using QGIS 1.8.0. Whenever I try to merge raster images I get the following error message:

Traceback (most recent call last):
File "C:\PROGRA~1\Quantum GIS Lisboa\bin\gdal_merge.py", line 36, in
import gdal
File "C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib\site-packages\gdal.py", line 2, in
from osgeo.gdal import deprecation_warn
File "C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib\site-packages\osgeo\__init__.py", line 21, in
_gdal = swig_import_helper()
File "C:\PROGRA~1\Quantum GIS Lisboa\apps\Python27\lib\site-packages\osgeo\__init__.py", line 17, in swig_import_helper
_mod = imp.load_module('_gdal', fp, pathname, description)
ImportError: DLL load failed: Die angegebene Prozedur wurde nicht gefunden.

Does anyone know what I sould do?


texts:
I have been trying out several solutions to try to work with GML 3.2 geometry served through WFS 2.0 data services.

There seems to be very little experience and software out there that manages to handle complex feature types and GML3.2 geometries.

I have been trying out 52 north's WPS package, which is supposed to fully support GML2.1 > GML3.2 and also should process data with plugins for GRASS, sextante etc.

However while trying to send a POST processing request, for example unsing one of my public WFS2.0 > GML3.2 geometries I get the following error:

  org.xml.sax.SAXException: Handler for gml:MultiSurface could not be found.

You can see this happening if you are sending a POST request to the following URL:
http://geoprocessing.demo.52north.org/latest-wps/WebProcessingService?

The GML3.2 data I am sending as input data:
http://inspire.biodiversity.ro/geoserver/ows?service=wfs&version=2.0.0&request=getfeature&typename=au:AdministrativeUnit&featureID=BC

or

http://inspire.biodiversity.ro/geoserver/ows?service=wfs&version=2.0.0&request=getfeature&typename=au:AdministrativeUnit&featureID=EPT

I am using the Open HTTP requester plugin from Firefox to submit this reqeust, sending the request as application/xml. This is the body of the request I am sending:

<Execute service="WPS" version="1.0.0"
xmlns="http://www.opengis.net/wps/1.0.0"
xmlns:ns="http://www.opengis.net/ows/1.1"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xmlns:xs="http://www.w3.org/2001/XMLSchema"
xmlns:base="http://inspire.ec.europa.eu/schemas/base/3.3"
xmlns:au="http://inspire.ec.europa.eu/schemas/au/4.0"
xmlns:br="http://inspire.ec.europa.eu/schemas/br/4.0"
xmlns:ps="http://inspire.ec.europa.eu/schemas/ps/4.0"
xmlns:gn="http://inspire.ec.europa.eu/schemas/gn/4.0"
xmlns:gco="http://www.isotc211.org/2005/gco"
xmlns:gmd="http://www.isotc211.org/2005/gmd"
xmlns:wfs="http://www.opengis.net/wfs"
xmlns:wps="http://www.opengis.net/wps/1.0.0"
xmlns:ows="http://www.opengis.net/ows/1.1"
xmlns:gml="http://www.opengis.net/gml/3.2"
xmlns:ogc="http://www.opengis.net/ogc"
xmlns:wcs="http://www.opengis.net/wcs/1.1.1"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:it.geosolutions="http://www.geo-solutions.it"
xsi:schemaLocation="
http://www.opengis.net/wps/1.0.0
http://schemas.opengis.net/wps/1.0.0/wpsAll.xsd
http://inspire.ec.europa.eu/schemas/au/4.0 http://inspire.ec.europa.eu/schemas/au/4.0/AdministrativeUnits.xsd
http://www.opengis.net/gml/3.2
http://inspire.biodiversity.ro:80/geoserver/schemas/gml/3.2.1/gml.xsd
http://www.opengis.net/wfs/2.0
http://inspire.biodiversity.ro:80/geoserver/schemas/wfs/2.0/wfs.xsd"
>  <ns:Identifier>org.n52.wps.server.algorithm.intersection.IntersectionAlgorithm</ns:Identifier>
  <DataInputs>
    <Input>
      <ns:Identifier>Polygon1</ns:Identifier>
      <Reference mimeType="text/xml; subtype=gml/3.2.1" xlink:href="http://inspire.biodiversity.ro/geoserver/ows?service=wfs&amp;version=2.0.0&amp;request=getfeature&amp;typename=au:AdministrativeUnit&amp;featureID=BC"/>
    </Input>
    <Input>
      <ns:Identifier>Polygon2</ns:Identifier>
      <Reference mimeType="text/xml; subtype=gml/3.2.1" xlink:href="http://inspire.biodiversity.ro/geoserver/ows?service=wfs&amp;version=2.0.0&amp;request=getfeature&amp;typename=au:AdministrativeUnit&amp;featureID=EPT"/>
    </Input>
  </DataInputs>
  <ResponseForm>
    <ResponseDocument>
      <Output mimeType="application/WFS">
        <ns2:Identifier xmlns:ns1="http://www.opengis.net/ows/1.1" xmlns:ns="http://www.opengis.net/wps/1.0.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ns2="http://www.opengis.net/ows/1.1">intersection_result</ns2:Identifier>
      </Output>
    </ResponseDocument>
  </ResponseForm>
</Execute>