texts: Intermittent 'Could Not Find Layer' and 'Bad Request' Errors on GeoServer

texts:
We are having some intermittent faults with GeoServer at the moment. We have a mapping application - http://explorer.geowessex.com - that makes WMS requests to our Geoserver to retrieve tiles. The layers are properly loaded and work via GeoServers built in 'Layer Preview' function, but sometimes in the application the tiles don't load. In the Firebug Net Panel you will see all the WMS Get Requests being fired when a layer is loaded and they say they return with no errors. In reality if you open one of these response in a new tab you get the following XML file.

<?xml version="1.0" encoding="UTF-8" standalone="no"?><!DOCTYPE ServiceExceptionReport SYSTEM "http://wms.geowessex.com:80/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd"> <ServiceExceptionReport version="1.1.1" ><ServiceException code="LayerNotDefined">
  Could not find layer 'workspacename':'layername'

Which is stupid because it definitely exists, I can see it load perfectly well in GeoServer layer preview. What's even more annoying is that sometimes it doesn't even return this, it just gives the error '400:Problem communicating with GeoServer'. On top of that I cannot reproduce it regularly, I can have the map up, turn on a layer and it wont work. Turn on another layer and it will, go back and turn the other layer back on and it will start working again.
I'm really banging my head against a brick wall here. GeoServers own documentation is not that good for troubleshooting and I must've done enough Google searches to bring down their servers in the last two days and am still no closer as to why this happens.

Anyone experienced this themselves?


texts: PCI Geomatica ATCOR error with RapidEye images

texts:
I'm doing atmospheric corrections to RapidEye images with ATCOR in PCI Geomatica 2014. However, there are serious salt-and-pepper type errors in the outcome bands 4 and 5. These errors only occur when haze removal is done as part of ATCOR. I would really need to do the haze removal as part of the atmospheric corection process.

Has anybody else come across to this error? Do you have any solutions to the problem?

Here is an example with RGB combination 5,4,3:


texts: Getting arcpy.AddMessage() to display message from script tool called from another script tool?

texts:
I am using a Python script tool within an ArcToolbox in ArcGIS 10.2.

The script calls another python script tool while running but does not display the messages from the called script when using the arcpy.AddMessage() method.

Example:

scribtbox (ArcToolbox)
two script tools

script1
script2

script1:

import arcpy

# output to the screen
arcpy.AddMessage("Doing routine 1...")
routine1()

# call script2
arcpy.script2_scriptbox

script2:

import arcpy

# output to the screen
arcpy.AddMessage("Doing routine 2...")
routine2()

The output model run dialog box in ArcGIS says:

> Doing routine 1...

and then finishes without saying the message from script2.

I'd like it to say:

> Doing routine 1...
> Doing routine 2...

and then finish.

Is this functionality not built in, am I doing something wrong, or is there some other issue at hand.


texts: How do I invert the color ramp in QGIS?  I am using RdYlGn, but I would rather use GnYlRd.  Seems really easy, but I can't seem to find where the invert is.  I am using 1.8.0 on Windows 7.


texts: How to get a string from a Double with a fixed number of decimals?

texts:
I have a featureclass (in a File Geodatabase), with a particular field of type Double. I want to copy it over to a text field, so that only 3 decimal places are shown.

For example: 34.78902345 becomes 34.789

I have gone through this question: How to convert a double field to a string?, but that does not solve my issue. I have even tried to format the double field, but even then the text field contains all the decimals.

How do I do this in ArcGIS 10.0?


texts: Assigning RouteIDs to polyline feature class in ArcGIS Linear referencing?

texts:
I have a Linear Referencing Route System in ArcGIS Geodatabase format (polylines with M values) and a separate polyline feature class of short segments (called TMC) that do not have RouteIds currently. The Linear Referencing Routes are directional. So, there will be two long lines on top each other in most locations in the Route layer where there are two-way routes, with RouteIds like 1001 and 1002 and these are digitized in opposite directions.

How can I assign the appropriate RouteIds to the TMC segment layer based on the digitized direction of the TMC segments?

I only need to assign RouteIds when the TMC is approximately parallel to an existing Route and within 100 ft of the Route? We need to ignore the TMCs which do not meet these criteria. So, if a TMC segment intersects a Route, but is perpendicular to it, it should not be assigned a RouteId.


texts: Unable to completely download OSM file to PostGIS?

texts:
I used this command to import the osm file to PostGIS:

osm2pgsql /Desktop/osm/mtl.osm -d demo_data -U blanco -P 5432 -S "/usr/local/share/osm2pgsql/default.style" --hstore

I am using osm2pgsql SVN version 0.84.0 (64bit id space)

Why am I unable to download planet_osm_nodes?

Using projection SRS 900913 (Spherical Mercator)
Setting up table: planet_osm_point
NOTICE:  table "planet_osm_point_tmp" does not exist, skipping
Setting up table: planet_osm_line
NOTICE:  table "planet_osm_line_tmp" does not exist, skipping
Setting up table: planet_osm_polygon
NOTICE:  table "planet_osm_polygon_tmp" does not exist, skipping
Setting up table: planet_osm_roads
NOTICE:  table "planet_osm_roads_tmp" does not exist, skipping
Using built-in tag processing pipeline
Allocating memory for dense node cache
Allocating dense node cache in one big chunk
Allocating memory for sparse node cache
Sharing dense sparse
Node-cache: cache=800MB, maxblocks=102400*8192, allocation method=3
Mid: Ram, scale=100

Reading in file: /Users/hadoop/Desktop/osm/mtl.osm
StartElement: Unknown element name: note
Unknown node type 3
EndElement: Unknown element name: note
StartElement: Unknown element name: meta
EndElement: Unknown element name: meta.
Processing: Node(1234k 112.2k/s) Way(248k 15.56k/s) Relation(6805 3402.50/s)        parse time: 29s

Node stats: total(1234539), max(4764987639) in 11s
Way stats: total(248884), max(483754469) in 16s
Relation stats: total(6805), max(7120865) in 2s
Committing transaction for planet_osm_point
Committing transaction for planet_osm_line
Committing transaction for planet_osm_polygon
Committing transaction for planet_osm_roads

Writing way (248k)
Committing transaction for planet_osm_point
Committing transaction for planet_osm_line
Committing transaction for planet_osm_polygon
Committing transaction for planet_osm_roads

Writing relation (6798)
Sorting data and creating indexes for planet_osm_point
Sorting data and creating indexes for planet_osm_line
Sorting data and creating indexes for planet_osm_polygon
Sorting data and creating indexes for planet_osm_roads
node cache: stored: 1234539(100.00%), storage efficiency: 52.39% (dense    blocks: 276, sparse nodes: 1036935), hit rate: 98.58%
Analyzing planet_osm_roads finished
Copying planet_osm_roads to cluster by geometry finished
Creating geometry index on  planet_osm_roads
Creating indexes on  planet_osm_roads finished
All indexes on  planet_osm_roads created  in 1s
Completed planet_osm_roads
Analyzing planet_osm_polygon finished
Analyzing planet_osm_line finished
Analyzing planet_osm_point finished
Copying planet_osm_polygon to cluster by geometry finished
Creating geometry index on  planet_osm_polygon
Copying planet_osm_line to cluster by geometry finished
Creating geometry index on  planet_osm_line
Copying planet_osm_point to cluster by geometry finished
Creating geometry index on  planet_osm_point
Creating indexes on  planet_osm_polygon finished
All indexes on  planet_osm_polygon created  in 5s
Completed planet_osm_polygon
Creating indexes on  planet_osm_line finished
All indexes on  planet_osm_line created  in 7s
Completed planet_osm_line
Creating indexes on  planet_osm_point finished
All indexes on  planet_osm_point created  in 8s
Completed planet_osm_point

Osm2pgsql took 54s overall


texts:
I have an Oracle database that I have been extracting the certain geometry tables into QGIS and then creating shapefile files for those tables.

Does anyone have any tips or tricks to automatic the process with python?

I also need to reproject the data before creating the shapefiles.


texts: How to get City, State, Numeration, etc from OSM roads by Geofabrik?

texts:
I'm really new in this GIS world. I downloaded a shp file from Geofabrik (http://download.geofabrik.de/osm/south-america/chile.shp.zip) and imported it to Postgresql using Postgis. If I see the road table, this is the info that I get:

http://grab.by/fHa4

Looking in internet, I understood that geom can give me some geographic info that I need (lat, long per example). But, I need to know more specific information about this polyline: City, State, numeration, etc. Is this info in geom column? Or in other part in this .shp? How could I get this information?


texts: OpenLayers 2.11 WMS layer in custom projection

texts:
I am using OpenLayers 2.11 from here http://cdnjs.com/libraries/openlayers
and trying to display WMS layer.
It is in 3301 projection so projection is defined as follows:

proj4.defs("EPSG:3301","+proj=lcc +lat_1=59.33333333333334 +lat_2=58 +lat_0=57.51755393055556 +lon_0=24 +x_0=500000 +y_0=6375000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs");
var projection = new OpenLayers.Projection('EPSG:3301');

and WMS layer and map like this:

 var bounds = new OpenLayers.Bounds(
            370753.1145, 6382922.7769, 739245.6000, 6624811.0577
    );
    var options = {
        maxExtent: bounds,
        maxResolution: 2601.4538070271437,
        projection:projection,
        units: 'm'
    };
    map = new OpenLayers.Map('map', options);
var base  = new OpenLayers.Layer.WMS("Maa/ameti aluskaart",
            "http://kaart.maaamet.ee/wms/alus?",{layers: 'MA-ALUS',
                version: '1.1.1' ,'srs':projection});
map.addLayer(base);

But there are no request at all with 2.11 version of OpenLayers. With OpenLayers 2.12 it is shown, but I am interesting in how to show WMS in older OpenLayers 2.11 version.


texts:
I am trying to create agent based model for simulating movement of occupants in different rooms of different floors of a building. I am using Agent Analyst tutorial downloaded from ESRI website and PhD dissertation of Elizabeth Groff as a guide.

Since the building has three different floors, I have created a 3D network dataset and created a toolboox which can find shortest route, closest facility and OD matrix from room centre to all the emergency exits.

I have created the following according to the two guides above.

csv file containing the activity space (i.e. the door of the room occupied, nearest exits of the floor, bends on the staircase and exit of the ground floor).
Point shapefile of path that the agents have to follow.
csv file containing the list of all neighbouring nodes.
Point shapefile of occupants. There are more than one occupants in few rooms. I intend to show the movement of these occupants as vector agents in Agent Analyst.
Point shapefile of all emergency exits. I also want it to be include in the model.

So I want to:

Show the movement of these agents from their present location to the nearest exit on ground floor.
The occupants are of different sex and age, thus I want the agents to move at different speed from one node to another and while they are moving downstairs. How can I do that?

Lastly, can anyone tell me if movement of agents can be shown in ArcScene because I tried to show the movement, it does not. I tried the same in ArcMap and it works. I haven't tried the same in ArcGlobe.

I am using ArcGIS 10 and Agent Analyst extension for ArcGIS.


texts: Batch convert ESRI Ascii Grid files to MapInfo TAB

texts:
I’m looking for a ‘quick’ (in terms of my time, not computer time) solution to converting a whole bunch of ESRI Ascii Grid files into MapInfo TAB files consisting of a grid of points with XYZ fields. I’ve looked at GDAL/OGR but as far as I can tell GDAL deals specifically with Raster/Grid data (i.e. the source data I’ve got) and OGR deals with vector data (my desired output). I can’t seem to find any command to go from grid to vector though.

I have the following software available: MapInfo Pro, Global Mapper v15, QGIS, GDAL

I'm aware that both Global Mapper and QGIS are scriptable but I've little to no experience in writing scripts for either of these so am avoiding that given the time it would take to a) learn what I'm doing and b) actually write the script.

I have considered implementing the conversion in MapInfo's scripting language, MapBasic, with which I am quite familiar, but it tends to be a bit slow and given the large number of files I have I'd rather find something a bit quicker. I also don't have a huge amount of time to spend on writing up a script for this so would rather not reinvent the wheel if something already exists.


texts: How to use itemChanged signal to update data?

texts:
I'm trying to understand the itemChanged signal to be able to update some information. An example of what I'm trying to achieve could be:

def setCoordText():
    label1.setText(str(label1.x()) + "," + str(label1.y()))

label1.itemChanged.connect(setCoordText)

so that when the label is moved to another location, it would display the new coordinates. What is wrong/missing here? Why nothing happens?


texts: How to unify the color of graduated classes (in symbology) across several maps?

texts:
I'm using ArcGIS 10.3

I have got 3 maps for the same parameter across several years.
I want the color of each class (e.g. 0.801-1.000) across all three maps to be the same to make comparison easily. This can be done manually by changing the color.

Is there any way to do it faster?


texts:
I am trying to set up Geoportal to harvest services hosted by Geoserver. I tried follow this instructions:

Open your Geoportal Server.
Go to the Administration tab Choose the 'Add' option
Choose 'Register resource on the network' and click Proceed
In the URL resource form, paste the OGC GetCapabilities URL (including all parameters)
Keep the options for allowing the resource to be found and for it to be synchronized at an interval
Choose the interval.
Select if you want to be notified of synchronization results.
Click Save.
In the Manage page, approve this newly registered resource.

Now I did the following in Geoserver:

Registered a ArcSDE shapefile in Geoserver that I can see in "Layer Preview" and I can get the common formats like OpenStreet map working for the layer.

But what is an OGC GetCapabilities URL??? Where do i declared and link it to the layer I registered?


texts: ArcGIS - How to copy and replace polygon (ONLY GEOMETRY)

texts: I have two layer (A, B). I need copy polygon from layer A and replace polygon in layer B, but only geometry. I want to preserve original attribute table from layer A.


texts:
I have two rasters in Geotiff, same CRS. I tried to merge them and I have this error code

ERROR 4: `C:/Anita_Tepoztln/Rasters/e142.tif' does not exist in the file system,
and is not recognised as a supported dataset name.

ERROR 4: `C:/Anita_Tepoztln/Rasters/e145.tif' does not exist in the file system,
and is not recognised as a supported dataset name.

Traceback (most recent call last):

File "C:\PROGRA~2\Quantum GIS Lisboa\bin\gdal_merge.py", line 509, in
sys.exit(main())
File "C:\PROGRA~2\Quantum GIS Lisboa\bin\gdal_merge.py", line 392, in main
ulx = file_infos[0].ulx
IndexError: list index out of range

I wonder what I am doing wrong.

I am able to load the rasters and see them, not merging them

I am using Qgis Lisboa standallone installation, on  Win7 64 16gb of RAM


texts: Assigning shortcut key to the Zoom to selection function is not working

texts:
I am trying to change the shortcut key cntrl + j (which is standard) to for example q, in the settings I see its changed but in the Attribute table where I want to use it to select a row and then zoom to selection its not working, cntrl + j is working.

How can I change this and use a different shortcut key?


texts:
ArcMap add-in is freezing on .Search() by spatial filter.
Here is my code, for ArcMap add-in 10.3 version. When I used it with Oracle databases, happens only on one instance, works perfectly with others. No error is thrown, it just waits there. Do you have any idea why this may be happening?

        IGeometry osakaGeometry = (IGeometry)point; //point - IPoint
        // Use the ITopologicalOperator interface to create a buffer.
        ITopologicalOperator topoOperator = (ITopologicalOperator)osakaGeometry;
        IGeometry buffer = topoOperator.Buffer(buffersize);

        List<IFeature> retVal = new List<IFeature>();
        ISpatialFilter spatialFilter = new SpatialFilterClass();
        spatialFilter.Geometry = buffer;//queryGeometry;
        spatialFilter.SpatialRel = reltype;//esriSpatialRelEnum.esriSpatialRelTouches;
        spatialFilter.SubFields = "*";
        spatialFilter.WhereClause = "1=1";

        foreach (IFeatureLayer fl in ToSelectLayers)
        {
            IFeatureCursor fcursor = fl.FeatureClass.Search(spatialFilter, false);

...
Not: I tried buffer size from 0.003 to 1. Apparently it is not the issue.


texts:
I am using Qgis2Web to create a webmap application. I have choose a basemap but it dosen't appears in the layerlist(collapsed layerlist).How could I edit the exported code (probably the qgis2web.js file ) in order to do this.

What i actually want to do is to add 2 or more basemaps in the final project and let the user choose from the layerlist.


texts: Symbolizing polygon fill color based on HTML hex color attribute value in ArcMap?

texts:
I have a set of polygons that I need to symbolize with a certain color. I have an column in the attribute table that contains the contains the html hex color value (#FF0000 Red) that I want the polygon to be colored. Is there a way to symbolize the color based on the attribute value?

I see that this might be possible using ArcPy along with a .lyr file, but I am not very experienced with ArcPy, and the last post isn't that promising, stating that the symbology property is read-only.
Using ArcPy to get layer symbology?

This post is making me think that it's not possible.
Changing color of symbol in UniqueValuesSymbology object using ArcPy?


texts: Writing the table layer_styles in a Geopackage

texts:
I made a pyQGIS script to backup the PostGIS tables of my team. I needed to save each table style in the Gpkg's layer_styles.
Here's an extract of how I've done but I would prefer to implement the pyQGIS one.
# Init QGIS
qgs = QgsApplication([], False)
qgs.setPrefixPath("C:/ProgGIS/QGIS314/apps/qgis", True)
qgs.initQgis()
#...

# PostGIS datasource
uri = QgsDataSourceUri()
uri.setConnection(myServeur, "5432", myBase, myUser, myPass)
uri.setDataSource('public', 'layer_styles', None)
# Importing the table layer_styles
couche_styles = QgsVectorLayer(uri.uri(False), 'layer_styles', "postgres")
couche_styles.setSubsetString(' "f_table_catalog" = \'myBase\' AND "f_table_schema" = \'mySchema\' AND "f_table_name" IN (\'table1\',\'table2\',\'table3\',\'table4\')')
QgsProject.instance().addMapLayer(couche_styles)
# Exporting the table layer_styles to the gpkg

#...
#exporting the other layers
#...

# Geopakagization of the gpkg's layer_styles (Empying id, f_table_catalog and f_table_schema)
connection_gpkg = sqlite3.connect(my_geopackage)
cursor_gpkg = connection_gpkg.cursor()
cursor_gpkg.execute("UPDATE layer_styles SET id = '', f_table_catalog = '', f_table_schema = ''")
connection_gpkg.commit()
cursor_gpkg.close()
connection_gpkg.close()

#...
# Quiting QGIS
qgs.exitQgis()


texts: Invalid expression in SelectLayerByAttribute?

texts:
I am coping with a shapefile and I want to select some rows if the values in one column are in a Python list. The following gives the code I write:

tpu_involved = list(set(dataframe['SmallTPU']))
selection_expression = '"SmallTPU" IN ({0})'.format(', '.join(map(str, tpu_involved)) or 'NULL')
print('---------------------------------------------------------------')
print('The invalid expression is: {}'.format(selection_expression))
print('---------------------------------------------------------------')
selected_layer = arcpy.SelectLayerByAttribute_management(in_layer_or_view='tpu_lyr',
                                                                 selection_type='NEW_SELECTION',
                                                                 where_clause=selection_expression)

The tpu_lyr is the layer which is loaded from a local shapefile. I have checked some previous questions like Selecting features by attribute if in Python list?, but it still does not work.


texts: Leaflet: How to move the layer-control menu?

texts:
This might be a silly question, but I could not find a documented way to achieve this.

I'd like to freely position the layer-control menu, probably to the top left near the default zoom-in/out button.

My layer control looks like this:

// Group layers as overlay pane
overlayPane = {
  "Endpoints" : endpointMarkerLayer,
  "Links" : linkLineLayer
};

// Add a layer control element to the map
layerControl = L.control.layers(null, overlayPane);
layerControl.addTo(map);

Where endpointMarkerLayer and linkLineLayer are layers containing markers and polylines respectively.

Is there an option to specify where the menu should appear?
Or alternatively, how can I get a reference to the DOM-objcet of the control-menu, such that I could assign it a custom class and override positioning in CSS?


texts: Optimize the locations of 2 additional restrooms on a trail (ArcGIS 10.2)

texts:
In the image below, the line is a trail and the point is an existing restroom. I would like to find the rough locations of 2 additional restrooms that minimize the maximum distance along trail from anywhere on the trail to the closest of the 3 restrooms.

The trail was GPS-collected in segments and connected with Edit Vertices in ArcMap. The dataset comprises 1 line shp for the trail and 1 point shp for the restroom. I have no idea how to proceed.

I have posted a simplified version of this question in stackoverflow.

Hints will be appreciated as well.
Thanks.


texts: How do I put label inside polygon at any scale in QGIS 1.8.0

texts:
I would like to put labels inside polygon anytime and any scale. I'm using "Labeling" from the toolbar. Although I tried all types of the placement from the advanced tab, some labels still fall outside of their respective polygon.

Please let me know how to put labels inside the polygon at any time and scale. Thank you.


texts: How to run an Overpass-API `around` query incrementally excluding already searched area?

texts:
I have to find some types of POIs around a given coordinates. Ideally, I want to find the nearest one of such POI around the given coordinates.
To meet my requirement I can search with either around
[out:json];
(
  node["public_transport"="platform"](around:500,48.2,16.3);
);
out body center qt;
>;
out skel qt;

or approximately with bbox
[out:json][bbox:48.19,16.29,48.21,16.31];
(
  node["public_transport"="platform"];
);
out body center qt;
>;
out skel qt;

My strategy to search quickly is to start search in a small area (short distance or smaller bbox), if I could not find anything then I will increase the search area. I search incrementally until something is found or a hard distance limit is reached.
I am not sure if Overpass-API implementation has done some performance tricks like cache of previous result.
Every time I increase search area I still search the previous smaller area. With bbox I can somehow split it and avoid already searched area, but I do not know I can achieve it with around.
How can I use around method to increase search area without covering old searched one?
update: to focus on search procedure itself. I will ignore the rate/resource limiting, latency, and so on, since I can eventually run my own instance instead of querying public one. I could perform queries in parallel or do some load balance among instance but it is just further optimisation options.


texts:
I am trying to compare two separate feature classes to identify differences between them (sort of a diff function). My basic workflow:

I extract the geometries using a SearchCursor
Save the geometries of the two feature classes as GeoJSON using a modified __geo_interface__ (got it from valveLondon return {'type': 'Polygon', 'coordinates': [[((pt.X, pt.Y) if pt else None) for pt in part] for part in self]}). This is to avoid the shared geometry object that ESRI uses with cursors and the inability to make deep copies (some discussions here on gis.stackexchange talk about it).
Check the geometries of the two feature classes based on a unique identifier. For example, compare the FC1 OID1 geometry with the FC2 OID1 geometry. To get the geometry as an ESRI object instance, call arcpy.AsShape() (modified to read polygons with holes (see point 2 above) with return cls(Array([map(lambda p: Point(*p) if p is not None else Point(), part) for part in coordinates])). The comparison is simply geom1.equals(geom2) as indicated in the Geometry Class.

I expect to find ~140 changes in geometries, but my script insists there are 430. I tried to check those GeoJSON representations and they are identical, yet the Geometry Class equals() refuses to say so.

An example is below:

>>> geom1geoJSON
{'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
>>> geom2geoJSON
{'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
>>> geom1 = arcpy.AsShape(geom1geoJSON)
>>> geom2 = arcpy.AsShape(geom2geoJSON)
>>> geom1.equals(geom2)
False
>>> geom2.equals(geom1)
False

The expected behavior here should be True (not False).

Does anyone have any suggestions before I move everything to ogr geometries? (I am hesitant as ogr.CreateGeometryFromGeoJSON() expects a string, and arcpy's __geo_interface__ returns a dictionary and I feel like I am adding extra complexity).

Found the following resources helpful, even though they do not answer the question:

arcpy.Geometry question here on gis.stackexchange.com which was linked above in my text.
Errors in arcpy's Polygon class from the arcgis.com forums (apparently there are a lot of precision errors in ArcGIS 10.0 which theoretically got fixed in 10.1 but I cannot verify that, in 10.0 SP5 you still get the error).


texts: Locating geometries within a polygon/area

texts:
I have a postgis database table of boundary data which is made up of LineString geometry's, what I would like to do is create a query that will return only the rows that contain a LineString geometry that is within (or at least one point of the LineString is within) a given polygon or area (this would typically be a bounding polygon of the current map view).

I've tried using ST_Intersects, ST_Contains and ST_Covers but to no success. What is the recommended method for doing this?


texts:
I am trying to create a map out of a file with millions of points. To make this tractable, I am using h3 to aggregate the points. After doing the aggregation using a combination of h3 and pandas in Python (adapted from here), I am converting the h3 indexes to polygon boundaries:
agg1=gdf.groupby(['h3', 'product', 'year']).agg(Total_Count=('product', 'sum')).groupby(level=0).apply(lambda x: 100*x / x.sum()).reset_index()

def add_geometry(row):
    points = h3.h3_to_geo_boundary(
      row['h3'], True)
    return Polygon(points)

agg1['geometry'] = agg1.apply(add_geometry, axis=1)

gdf_op = gpd.GeoDataFrame(agg1, crs='EPSG:4326')
output_filename = 'gridcounts.gpkg'
gdf_op.to_file(driver='GPKG', filename=output_filename)

However, when I open the file in QGIS, some hexagons are streaking across the map:

I have a feeling it is because h3 is defined on a globe and as a result the polygons towards the edges translate to coordinate bounds which are at the opposite ends of the projected map. Does anyone have a solution to this apart from deleting the problem hexagons?


texts: I am unable to write the closing square bracket in QGIS Python console or editor, unless automatic parentheses insertion is enabled. Parenthesis and curly brackets work fine. Keyboard layout is Spanish. Happens only in QGIS. Tried two different computers. Is this a bug or what is going on? Basically, I can't type " ] " in QGIS Python


texts: ArcGIS map service and image service in Cesium

texts:
Im trying to load an ArcGIS map service in Cesium with WMS protocol but I get this exception:

Error: error occurred in "v": Failed to obtain image tile X : 0 Y : 0 Level : 0.

When I publish the map service without WMS protocol it works fine (using ArcGisMapServerImageryProvider). I get the same exception when trying to show an image service.

I tried use both WebMapServiceImageryProvider, that is suppose to support WMS protocol, and ArcGisMapServerImageryProvider but that did not work as well.

For the map service I use this url: http://server/arcgis/services/MapName/MapServer/WMSServer. For the image service I use this url: http://server/arcgis/services/MapName/ImageServer/WMSServer

Could anyone help me show both map service published in WMS protocol and image service?


texts: Generating an invalid vector layer in QGIS

texts:
I need an invalid layer to test a plugin in QGIS. In the plugin window, I select a layer from combobox. So, the invalid layer should be in Layers panel. I google for a while but I cannot find a way to generate an invalid vector layer.
How can I generate one so that layer.isValid() returns False?


texts: Performing aerial triangulation on scanned aerial photos without having surveyed ground control points (GCPs)?

texts: I have a set of 1200 historical aerial photos covering North Carolina outer banks area from 1974, 1984 and 1992. I need to georeference them and produce seamless color balanced dataset that will be utilized for compilation of NC outer banks atlas. Except for scale information, flight number and date on photo paper copies, no other information regarding camera positions or ground control points is available. No matter how well I pick control points for image-to-image rectification (I use ArcGIS base map for that), the photos still don't align properly when I display them together in ArcMap. I learnt that analytical aerial triangulation may significantly improve the results, but it seems like it needs GCPs that I don't possess. So, my question is how can I triangulate my photos (or fix my alignment problem in any other way) without GCPs and preferable using ArcGIS, Agisoft Photoscan or other open source software.


texts: Is it possible to change the field length of a shapefile to over 254?

texts:
The maximum field length for a shapefile is 254 (https://en.wikipedia.org/wiki/Shapefile#Limitations) and because of this I'm having problems with the input attributes, which have a size greater than 254.

What alternatives are there for having longer field lengths?


texts: QGIS layout template with embedded pictures

texts: I have designed different layout templates with a couple of logos (pictures) in QGis 3.12.0, in order to share these with some project partners. The maps can be loaded via PostgreSQL database connection. The problem is that the logo-files (picture.jpeg) are  loeded from my local computer and get lost when the project partners are loading the qps-file or there own local Qgis-project. Is there a way to embedd the pictures in the layout file?


texts: Prompting for user input during ArcGIS Python script tool execution?

texts:
I want to prompt for user input after all parameters have been specified for the script tool. I've seen many forum posts and stackexchange questions about this, and they all conclude that it's unnecessary, but I think my case is an exception. I want to give the user an opportunity to confirm data, such as how many sets of GPS coordinates there are. They will be yes/no questions ("is this data obviously incorrect?"), and the script tool parameter input form doesn't seem like it'd work well for this.

I know Tkinter's event handler causes problems with ArcGIS, and other framworks have similar problems. One workaround I found is to have your script call another script that runs as a separate process, and display the GUI from the second process, but that seems a little too hackish.

Apparently wxPython does work, but it's not included with python. I need to be able to distribute these scripts without any extra requirements other than ArcGIS, so that wouldn't work for me.

To clarify, the most important use would be to prompt the user when my software validation checks can't validate the data. For example, I have a list of IDs that I want to validate, but the format isn't consistent. I want to be able to prompt the user with a yes/no dialog if I can't validate in software (eg. the ID doesn't match one of the common patterns). It's possible that new ID formats will be created without my knowledge.

Another example: I have a file that contains, among other things, a table of start & stop times. I want to validate each row of the table individually, and cramming all of that data, plus the other checks, into the warning window would get ugly fast.


texts: Calculate bearing between two decimal GPS coordinates (Arduino / C++)

texts:
I am trying to calculate the bearing just like in this example: Calculate bearing between two decimal GPS coordinates
, but my environment is Arduino.

I got the code translated (please see below), but I am having issues with the floating points. I need a higher resolution/precision.

Using these points:

double r1 = 9.935483816921979;
double r2 = -84.05334153306967;
double g1 = 9.936029735119176;
double g2 = -84.0532617256946;

I get these results:

8.194077303640654 // In the js code from the example
9.229920 // from my Arduino code

How can I get more precision?

// Here is my Arduino Code
double getBearing(double lat1,double lng1,double lat2,double lng2) {
  lat1 = radians(lat1);
  lng1 = radians(lng1);
  lat2 = radians(lat2);
  lng2 = radians(lng2);

  double dLng = lng2 - lng1;
  double dPhi = log(tan(lat2 / 2.0 + PI / 4.0) / tan(lat1 / 2.0 + PI / 4.0));

  if (abs(dLng) > PI){
    if (dLng > 0.0)
     dLng = -(2.0 * PI - dLng);
  else
    dLng = (2.0 * PI + dLng);
  }

  return fmod((degrees(atan2(dLng, dPhi)) + 360.0), 360.0);
}


texts: How do I get UTM projected areas from a shapefile created in QGIS?

texts:
I create a polygon (or point or line) shapefile in QGIS and it is automatically saved as WGS84 (as I understand, correct me if I'm wrong).  That is fine by me despite being a bit inconvenient.  Then I need the area.  I know that QGIS can generate areas and update the attribute table but it is unclear what the units are and WGS84 units don't readily translate to working area units such as acres.  I messed with this for a while and decided to work in Grass since I didn't know what sort of units I was getting and saving the layer in different CRS designations didn't result in different values in the attribute table.

I try to import the QGIS generated shapefile to Grass so I can reproject it to UTM and I get an error about illegal values.  The Grass environment was created using the QGIS generate file so it seems to me that QGIS is creating incompatible files when I create polygon shapefile.

Does anyone have any insight about this and how I can get legal values from which I can project to my desired projection?


texts: Making 3d block from flat 2d polygon in ArcScene?

texts:
I would like to make a mountain from flat polygon (this polygon is the base contour of the mountain). The question is how to make it keeping 35 angle of the slope ? The top of the mountain should be flat (no cone). At the end I would like to make calculation of the volume.

I'm new in ArcScene.


texts: Using Python to return elevation at arbitrary XY, from contour vector file

texts:
I want to write a Python3 script to calculate the elevation difference between a 3D surface (supplied as a DXF of contours) and a long list of XYZ coordinates (supplied as a CSV).  The script should have a minimal requirement for additional modules or installed applications.
The logical steps are easy:

Import one or more suitable Python modules such as GDAL or GeoPandas
Write a function elevation(X,Y) to return an elevation at a given XY coordinate.  This might have a prerequisite such as generating a vector TIN or a raster DEM.  For speed, generalizing the contour layer might also be a useful pre-process.
Open the contour DXF and perhaps post-process it to provide the "feed" needed for the elevation function
Loop through the data from the CSV, outputting the elevation difference for each row.

I'm having problems finding a way to do this job.  Research keeps leading to graphical GIS environments, which should not be necessary.


texts: Why are heatmaps on qgis 2.8.1 produced in the wrong CRS?

texts:
I noticed that the heatmap calculation on the new upgrade version of qgis is given the raster map without any coordinate system.

Have you got the same problem?


texts: Is there limit to number of elements in layout of ArcGIS Pro project manipulated by ArcPy?

texts:
Has anyone encountered, or aware of, a limit to the number of elements in an ArcGIS Pro layout that is manipulated by ArcPy?

The reason I ask is that I have created an ArcGIS Pro 1.4.1 project that has one layout and six map frames that I have been happily running a large ArcPy script against that finishes by exporting to PDF and then saving a copy of the project.

As part of the code I incorporate writing a table that incorporates cloning hundreds of text and graphics elements, and this is based on TextElement example 2 from the ArcGIS Pro documentation which:

  constructs a graphic table based on data values from a table in a map.
  The layout was authored with a vertical line named vertLine, a
  horizontal line named horzLine, and a text element named cellText.
  Each of the elements were authored with the appropriate symbology and
  text properties. The element's anchors were also set to the upper left
  position, and the text element's vertical and horizontal justification
  were set to upper left.

That works fine but when I add more code to implement a workaround that I am developing for measured grids (because they are not supported until ArcGIS Pro 2.x) that also creates about a hundred text elements, this new set of text elements do not appear in the exported PDF or the project saved as a copy, but using print statements in the code I can see that they are present in the original manipulated project.

As soon as I comment out the TextElement example 2 code, the measured grid text elements start to appear in both the exported PDF and the project saved as a copy.  This leads me to theorize that I have encountered a limit to how many layout elements ArcGIS Pro can support when one of its projects is manipulated by ArcGIS Pro.

I have not yet included a code snippet here, in case someone already has one that they can offer, but I believe I will be able to run up a code snippet that illustrates this.


texts: ArcGIS Runtime SDK for Qt - Drawing Distance Line

texts:
I'm just starting with ArcGIS Runtime SDK for Qt and having gone through many of the examples, I'm looking for advice as to how I would draw a line representing the distance between two points.

My understand thus far is that I can do this in either a graphical layer or in a graphics overlay but which is the preferred option?

I have managed to capture the start and end coordinate in WGS-84 but in terms of doing the actual drawing, I'm a little lost at present


texts:
I'm trying to create a custom layer, nothing crazy just press a button and have something show up on the display and also have the layer display in the TOC.

Well, the layer shows in the TOC but nothing shows on the display and the resulting error in my Debug tells me

  A first chance exception of type 'System.Runtime.InteropService.ComException' occured in MapControlApplication1.exe

I'm going with fairly straight forward code.

I've created a BaseCommand that ties into the toolbar.  When I press the button on the toolbar it successfully goes to Onclick.  Here is the meaningful code in the OnClick

CustomLayer2  layer = new CustomLayer2();
layer.Name = "Custom Layer";
mapControl1.AddLayer(layer , 0);

Here is the code for CustomLayer2

public class CustLayer2: BaseCustomLayer
{
    private IPolyline m_Polyline = null;

    public override void Draw(esriDrawPhase drawPhase, IDisplay Display, ITrackCancel trackCancel)
    {
        if (drawPhase != esriDrawPhase.esriDPGeography)
            return;
        if (m_Polyline == null)
        {
            IProjectedCoordinateSystem pcs = Display.DisplayTransformation.SpatialReference as IProjectedCoordinateSystem;
            if (pcs != null)
            {
                m_Polyline = MakePolyline(pcs.GeographicCoordinateSystem);
                m_Polyline.Project(pcs);
            }
            else
            {
                m_Polyline = MakePolyline(new UnknownCoordinateSystemClass());
            }
        }
        Display.DrawPolyline((IGeometry)m_Polyline);
    }

    public static IPolyline MakePolyline(ISpatialReference sr)
    {
        IPointCollection pc = new PolylineClass();

        object missing = Type.Missing; // (not needed with VB.NET)
        pc.AddPoint(MakePoint(-98.0, 29.0, sr), ref missing, ref missing);
        pc.AddPoint(MakePoint(-97.0, 28.0, sr), ref missing, ref missing);
        pc.AddPoint(MakePoint(-96.0, 27.0, sr), ref missing, ref missing);
        pc.AddPoint(MakePoint(-95.0, 26.0, sr), ref missing, ref missing);

        IPolyline polyline = (IPolyline)pc;
        polyline.SpatialReference = sr;
        return polyline;
    }
    public static IPoint MakePoint(double x, double y, ISpatialReference sr)
    {
        IPoint p = new PointClass();
        p.PutCoords(x, y);
        p.SpatialReference = sr;
        return p;
    }
}

Any ideas on why this would error out?  Note that a new layer is shown in the TOC but obviously the Draw never finishes.  The error is occuring on the Display.DrawPolyline((IGeometry)m_Polyline);


texts: Moved labels in QGIS 1.8 disappear on saving vector file

texts:
From other posts I worked out how to activate the "move label" icon and was successfully able to move a number of labels into the positions I wanted. However when I untoggle the "edit" function for the vector file concerned the labels that I have just moved all disappear - leaving me with only the unmoved labels.

I cannot find anyway to bring them back. How do I stop the moved labels from disappearing?


texts: Plot shapefile on top of raster using plot and imshow from matplotlib

texts:
I need to plot a shapefile on top of a raster. I would like to use the python package rasterio and some combination of fiona or geopandas to load the shapefile.

I found an example in the rasterio documentation but it doesn't provide code for the actual plotting.

I have tried loading the raster with rasterio, loading the shapefile using fiona, calling imshow on the raster, and calling plot on the shapefile. However, the shapefile does not appear on the plot.


texts: EsriLoader in Angular: basemap without border and streets

texts:
I'm using the Esri maps inside my angular 9 application by implementing the Esri loader like described in the tutorial. The data requires that there are no borders displayed (operas from 17xx where correct borders are hard to obtain).
Is there a basemap without borders, streets and names of countries? Just a gray shape to distinguish between ocean and land would be enough. City names would also be okay.
Because I wasn't able to find one, I played around with the ArcGIS online tool and found out that using the National Geographic and setting the transparency of National Geographic Style to 100% would fulfill my needs. How can I set the transparency of this layer in my application to 100%?


texts: QGIS- Try writing plugin - neither 'make' nor 'pyrcc4' is known command

texts:
I want to write a python plugin for QGIS 3.6.

Following the Guide I've installed the mentioned software Qt Creator and the  OSGeo4W . There is written, that in the OSGeo4W Shell It would be possible then to call the pyrcc4 command.

Later in the tutorial done by "make pyrcc4 ..."

However, nothing works on my machine.
No matter what I want to call, neither 'make' (like in Tutorial) nor 'pyrcc4' (like in some helps like this: Making 'make'-command work for compiling QGIS-plugins?) is a known command using the OSGeo4W Shell.

I've also tried this:

pyrcc4 command not recognized as internal or external command?

What am I doing wrong?


texts: Clarification of Logical Operators in a PostGIS Database

texts:
The following query returns this result:

select count(*) from a , b where a.geom && b.geom;

     count
    -------
     66514

On a postgis 2 database.

What does the logical operator AND (&&) between two geometry columns return?

What about the other operators?

Edit: After a quick test it seems it returns all the rows where the bounding boxes of both the geometries overlap with each other


texts:
I want to calculate isochrones representing the accessibility of certain locations in a city. So I have set up the pgrouting enabled database, imported the OSM-Roads and created a network. I have also calculated the "drivetime" of the road segments (time for walking, biking) and created new columns in the network table for this.

Now I have my point shapefile with the locations, which are subject to the accessability research (see screenshot).

My questions are:

How do I connect these location points from my point shapefile, which i imported into the database to the road network (I assume this is necessary in order to create the isochrones)?
Then how do I create 5min, 10min and 15min isochrones (for walking and cycling) for each of my 15 points (locations)?

I work on this for weeks now and also tried GRASS where I was at least able to connect my point shape points to the network but the outcome of v.net.iso is not what i am looking for.

A similar problem was discussed in How to connect points (bus stops), which don't lie on the lines (LINESTRING), to a network? but no solution was found. Also this popular example of finnish airports does not help me because I cannot see how the airports were "snapped" to the network from the explanation.

Could anybody help me with this?


texts: Working with feature classes from multiple geodatabases in ArcPy?

texts:
I have produced a Python script which creates multiple buffers of specific (listed) feature classes within a geodatabase. This is achieved by running each feature class through a for loop.

I would like to know if/how it would be possible to list feature classes from more than one geodatabase?

Currently I am using env.workspace to specify the single geodatabase. Listing the feature classes and then running them through the 'for loop'. A possible problem might be that many of the feature classes in the various geodatabases have the same name.

import arcpy
import os.path
import tempfile
from arcpy import env

# sets the workspace to my gdb
arcpy.env.workspace = "C:/GIS Home/Mining Features (MATTHEW).gdb"

# Output workspace and temp location is defined
output_workspace = "C:/GIS Home/project_1"
temp_workspace = tempfile.mkdtemp()

# list of features to buffer in gdb
fcList = ["Adit_Metalliferous", "Shaft_Metalliferous", "Quarry_Metalliferous", "Dump_Metalliferous", "Engine_House", "Goffan_Working", "Lode_Surface_Outcrop", "Lode_Undefined_Elevation", "Mine_Level", "Mine_Raise", "Portal", "Reported_Subsidence", "Secured_Feature", "Tunnel"]

#create empty lists for output
l_250 = []
l_500 = []
l_750 = []

for featureClass in fcList:
    b250 = os.path.join(temp_workspace, featureClass + "_b250.shp")  # make temp output file
    b500 = os.path.join(temp_workspace, featureClass + "_b500.shp")  # make temp output file
    b750 = os.path.join(temp_workspace, featureClass + "_b750.shp")  # make temp output file
    arcpy.Buffer_analysis(featureClass, b250, 250, "", "", "ALL")   # carry out buffer
    arcpy.Buffer_analysis(featureClass, b500, 500, "", "", "ALL")   # carry out buffer
    arcpy.Buffer_analysis(featureClass, b750, 750, "", "", "ALL")   # carry out buffer
    l_250.append(b250)  # append the buffer to the l_250 list
    l_500.append(b500)  # append the buffer to the l_500 list
    l_750.append(b750)  # append the buffer to the l_750 list

merge_250 = os.path.join(temp_workspace, "250_merge.shp")
arcpy.Merge_management(l_250, merge_250)
merge_500 = os.path.join(temp_workspace, "500_merge.shp")
arcpy.Merge_management(l_500, merge_500)
merge_750 = os.path.join(temp_workspace, "750_merge.shp")
arcpy.Merge_management(l_750, merge_750)

dissolve_250 = os.path.join(output_workspace, "250_dissolve.shp") # create 250_dissolve.shp for final output
arcpy.Dissolve_management(merge_250, dissolve_250)
dissolve_500 = os.path.join(output_workspace, "500_dissolve.shp") # create 500_dissolve.shp for final output
arcpy.Dissolve_management(merge_500, dissolve_500)
dissolve_750 = os.path.join(output_workspace, "750_dissolve.shp") # create 750_dissolve.shp for final output
arcpy.Dissolve_management(merge_750, dissolve_750)


texts: Problems with conversion of multilinestring to linestring

texts:
I'm trying to convert a column of datatype multilinestring to linestring using sql:

ALTER TABLE streets_nordeste
ALTER COLUMN geom TYPE geometry(linestring,3857) USING ST_GeometryN(geom, 1);

But the result is:

ERROR: Geometry SRID (0) does not match column SRID (3857)
SQL state: 22023


texts: (DSAS) Baseline won't Transect Appended Shoreline

texts:
I am trying to Transect my "Baseline" with my appended shoreline using DSAS in ArcGIS version 10.3. But this message window keep popping up every time I transect ( see attached first photo).

Its telling me to check out the DSAS log.txt so when I check it out this is what it says (see the 2nd attached photo) :

The geometry is legal but certain of its properties cannot be computed.

I think this is the problem but I can't seem to identify what the problem is specifically or what really went wrong. I thought maybe I did something wrong with my appended shoreline so I edited the  values, changed the Date configuration, everything I could possibly think of. But the message window I was talking about earlier just keep popping up. What do you think should I do?


texts: How to merge many shape with the same name in different folders?

texts:
I've many (hundreds) shape (3D points and 3D polygons) that have the same name, saved in different folders; and I've tried to merge them using batch with Merge tool in ArcGIS... but inexplicably ArcGis say that cannot process file with the same name! I'Now it's clear that I've too many files to rename them.

I've then used QGis and the process works, but lose the elevation information (Z) that I need.

I've seen other analogue questions but nothing directly usable. Some suggestions for a simply solution?


texts: Remove inner rings from QgsGeometry polygon in QGIS Python

texts:
I am a QGIS python user.
I want to iterate over the rings of a polygon to remove holes in the polygon, so:

I need to iterate over the rings in a QgsGeometry polygon
I need to detect the main ring of the polygon
I need to remove every ring, which is not the main ring

How can I do that in python QGIS?


texts: How to save a QGIS 2.0 project template?

texts:
When I attended the QGIS in Action workshop at FOSS4G13 Tim Sutton mentioned there was a new option to save a project as a template which when you opened up a new QGIS project you could pick to choose the saved template.

In my QGIS 2.0.1 Windows from the Project menu option there is a "New from Template" but how can I save a project as a template?

When I choose Save As... it is just a standard project file.

This may be linked with Print Composer but I have not found any guidance on using this option.

Thanks


texts: How to remove popup window from KML layer in ArcGIS API for JavaScript?

texts:
I am using the ArcGIS API for JavaScript for creating a map. I have added some CSV and KML layers to the map.

Now the KML is colored base layer(different colored for country) on the map and CSV is some green dot layer on map but top on the KML layer.When I click on the green dot(CSV) the green dot is selected and a popup window is shown on the for green dot(CSV layer).But on clicking the green dot  the base KML layer is also selected and also give popup window for the KML layer on goto the arrow sing on the header of popup.

Here is a link to my demo. You can check here click on the green dot and click on the arrow sing on the right of popup  black header.

Can anyone suggest how to remove the popup from the KML layer?


texts:
I'm using QGIS for quite a while now but I'm having some problems in measuring the volume of a depression zoned North of a waste dump I'm studying.

For the analysis I'm using a recent 1mx1m DTM of the area:

As I said the main element of the study is an estimation of the volume of the depression North of the dump. Here's the hillside of the area:

where I highlighted the boundaries of the depression.

As tool for the estimation I used the SAGA GIS tool Raster Volume in Qgis.

As Input GRID information I used the 1mx1m DTM but after running it on the log message panel here's the result:

Usage: saga_cmd grid_calculus 2 [-GRID <str>] [-METHOD <str>] [-LEVEL <str>]
            -GRID:<str>     Grid
            Grid (input)
            -METHOD:<str>   Method
            Choice
            Available Choices:
            [0] Count Only Above Base Level
            [1] Count Only Below Base Level
            [2] Subtract Volumes Below Base Level
            [3] Add Volumes Below Base Level
            Default: 0
            -LEVEL:<str>    Base Level
            Floating point
            Default: 0.000000*

I used the "Count only above base level" method considering the lowest point of the DTM as base level, but as you can see the result is a 0 value volume.

Am I mistaking anything? Here is the .asc DTM.


texts: GeoGratis API, 'OR'ing categories in a single request

texts:
I am a developer trying to work with the GeoGratis web API, provided here: http://geogratis.gc.ca/api/en

I want to make a request like this:

http://geogratis.gc.ca/api/en/nrcan-rncan/ess-sst/-/OR/(urn:iso:type)map/(urn:gc:subject)climate

but no response is given. This request does work:

http://geogratis.gc.ca/api/en/nrcan-rncan/ess-sst/-/AND/(urn:iso:type)map/(urn:gc:subject)climate

The difference between the two requests is the use of 'OR' in the first and 'AND' in the second.

I would like to be able to make one request that returns results from multiple groups, currently it looks like I have to make multiple requests and collate the results manually.


texts: Does somebody know where you can find the relative orbit number of a Sentinel 1 image? I know it should be somewhere in the metadata (the xml-file), but I can't find it.


texts: Comparison of clustering methods on geodemographic classification

texts:
Most geodemographic classification systems using K-means method for clustering and creating segments. For example, ESRI Tapestry(Business Analyst) is most using geodemographic system lately.

What about other clustering methods like CAMEO, ACORN, or Open geodemographics? Is there any way to compare clustering methods for geodemographics?


texts: OpenWind Unable to open polygon shp files from QGIS

texts:
I have a set of interlinked lines shp files which i imported to QGIS and add a fixed buffered zone via vector processing function on QGIS 3.10.
I saved the scratch layer and toggled edit off on polygon layer and then try to open this save layer in Openwind but got this error message: "Record found which is a different type (0) to that in main header (5) ERROR: Unexpected end of text."
Can i please ask for some help to see what have i done wrong and how can i correct this issue? This does not appear to be an issue with i deal with point layers so far (also used QGIS to add buffer)
"


texts: How to dynamically change MapServer's GRID MININTERVAL?

texts:
I started looking into using something similar to this example of the GRID layer definition in a MapServer mapfile that will draw lat/lon line grid at least every 5 degrees:
LAYER
  NAME "grid"
  METADATA
    "DESCRIPTION" "Grid"
  END
  TYPE LINE
  STATUS ON
  CLASS
    NAME "Graticule"
    COLOR 0 0 0
    LABEL
      COLOR  255 0 0
      FONT "sans"
      TYPE truetype
      SIZE 8
      POSITION AUTO
      PARTIALS FALSE
      BUFFER 2
      OUTLINECOLOR 255 255 255
    END
  END
  PROJECTION
    "init=epsg:4326"
  END
  GRID
    LABELFORMAT '%g°'
    MININTERVAL 5
  END
END # Layer

However, I don't want the MININTERVAL value to be hardcoded and would like to be able to vary that value by entering a different value in the URL call from a web page using the map_object notation.  But, I've been unsuccessful.  I have tried this in the URL call to MapServer:
&map_layer[grid].grid[0]=MININTERVAL+10 #changing to 10 degrees instead of 5
using the second example at this link: https://mapserver.org/cgi/controls.html#cgi-mapfile-change-parameters.  But, I either get a parsing error or an Internal Server error.  I've tried turning debug level to 5 and looking at logs but they don't offer any help.  Perhaps I'm just not using the correct map_object notation, but I'm not sure what else to use given the examples in the link.  I also attempted to look at the MapServer code to no avail.
I'm using MapServer v 7.2.2 but it seems like the GRID layer has been available for a while now.


texts: Open Basemap in different Projection (Transverse Mercator) - QGIS

texts:
I have data in a transverse Mercator projection which is a different projection from my base map (Spherical Mercator / Web Mercator / pseudo-Mercator). I really would like to create a map in a TM projection (using open/free basemap) but I couldn't find any open basemap with Transverse Mercator projection.

How can I create a map in TM projection where I include a base map in QGIS?
Is there any base map that provides a different projection than Web Mercator?

Note: I'm using QGIS. In Arcmap this is not an issue since they provide TM base maps. If you zoom in (greater than 1:6000000) it ok but smaller than that you'll find your basemap pixelized.


texts: Bad PostGIS TIGER Geocoding Performance In Specific States But Not Others

texts:
Short version: What should I do when I get great geocoding performance for all but one installed state?

Long version: I (not a DB expert) am running the PostGIS tiger geocoder on my machine, not a fancy server or anything. For addresses in two of the states, IL and NY, performance is awesome. Getting 20-40 ms per address (using EXPLAIN ANALYZE) when I just test the plain geocode function on a string. The default batch script works fine for these states as well.

For CA (using EXPLAIN ANALYZE) I get 1-2 minutes. At first I thought it was something in the batch script or a problem in the (larger) source address table, but the relative problem persists when geocoding a single hand entered string.

Here are the results for EXPLAIN on the basic geocode function.

explain analyze select geocode('ILLINOIS ADDRESS STRING')

Result  (cost=0.00..5.25 rows=1000 width=32) (actual time=23.943..23.944 rows=1 loops=1)
Planning time: 0.025 ms
Execution time: 23.954 ms

explain analyze select geocode('[NEW YORK ADDRESS STRING]')

Result  (cost=0.00..5.25 rows=1000 width=32) (actual time=39.697..39.698 rows=1 loops=1)
Planning time: 0.024 ms
Execution time: 39.711 ms

explain analyze select geocode('CA ADDRESS STRING')

Result  (cost=0.00..5.25 rows=1000 width=32) (actual time=132202.375..132202.375 rows=1 loops=1)
Planning time: 0.025 ms
Execution time: 132202.388 ms

This seems to indicate that there's some problem with the tiger data for California maybe? I really have no idea where to start to address this. I installed and configured this DB a month or two ago and I think I remember successfully testing some California addresses at this time. That suggests something has happened since I installed it, or that I'm misremembering a detail.


texts: How to get JSON response from submitJob operation?

texts:
I'm trying to get the JSON response of the submitJob operation launched by a restful call: http:///submitJob.
I'm using Java7.

The response I'm getting is an HTML response, using the code below:

URL url = new URL("http://<gp-task-url>/submitJob");

HttpURLConnection connection = (HttpURLConnection) url
        .openConnection();
connection.setRequestProperty("Content-Type", "application/json");
connection.setRequestMethod("GET");

InputStream is = connection.getInputStream();

Do you have any idea on how to retrieve a JSON response?


texts: How to rotate u and v components of wind in rotated gribs?

texts:
I have a grib message of a rotated hirlam product.

The grib message states, that the u and v component contained are relative to the rotated grid.

How can I unrotate them?


texts: PostGIS 2.1 error install with shared object (rtpostgis-2.1.so") Linux Mint

texts:
I've Installed Linux Mint 16 mate and install PostGIS with PostGreSQL (I've tried many ways in the Internet from original source to some suggestion) but all is too error: "
ERROR:  could not load library "/usr/lib/postgresql/9.3/lib/rtpostgis-2.1.so": libxerces-c-3.1.so: cannot open shared object file: No such file or directory"

Please help me for fix this problem because I can't understand how to fix this. Please, note that in "/usr/lib/postgresql/9.3/lib" has file rtpostgis-2.1.so.

Thanks.


texts: Is there a way in arcpy to convert map units from feet to decimal degrees?

texts:
This is my problem:
I have an arcpy script that finds the XMax, XMin, YMax, and YMin of a polygon. It shows the results in feet, but I want them to be in decimal degrees.

In ArcMap I can change the data frame's display units to decimal degrees. But every time I run my script it shows the numbers in feet again.

Is there a solution for this in arcpy?


texts: Listing Tables for several MXD documents using ArcObjects?

texts:
I want to do a "Extending ArcObjects: Windows Application" app:

I have a code snippet that lists all the layers from the MXD files stored in a directory:

string[] fileEntries = Directory.GetFiles("C://path");
foreach (string fileName in fileEntries) {
            IMapDocument mapDoc = new MapDocument();
            mapDoc.Open(fileName);
            IMap map = mapDoc.Map[0];
            for (int i = 0; i < map.LayerCount; i++)
            {
                ILayer l = map.Layer[i];
            }

That part is working perfectly.

Now I need to list all the tables, but I don't know how to create the IMxDocument object I need for that. It looks like that with that method I can only get the Tables from the current document in an Add-In.

How can I get the Tables of the MXD files stored in a directory using ArcObjects?

I already know how to do this with ArcPy, but I really need to be able to do the same with ArcObjects.


texts: Identifying Coordinate System of Shapefile when Unknown?

texts:
I have a Shapefile but its coordinate system is Unknown, and there is no *.prj file.

How can I identify it now?


texts: How to fix "Access to database denied" when accessing PostgreSQL?

texts:
I got a tiny problem regarding connection to my server. What should we do when this message appears in front of you?

I've done several changes in the pg_hba.conf entry but yet to succeed. Any careless mistake from my side?


texts: MapProxy doesn't start with EPSG:28992 while it starts and works fine with EPSG:4326

texts:
We are moving from Debian to CentOS and in that process we discover that MapProxy installation is turning out tougher in CentOS than it was in Debian.
One of the dependencies 'libproj' is not available as a CentOS package so I installed 'proj' which seems similar.

But when I try to start the MapProxy Instance it gives the following error stack:

[xyz@localhost spinmapproxy]$ sudo mapproxy-util serve-develop mapproxy.yaml
[sudo] password for xyz:
[2019-03-25 10:22:31,154] mapproxy.config - INFO - reading: /opt/spinmapproxy/mapproxy.yaml
Traceback (most recent call last):
  File "/bin/mapproxy-util", line 11, in <module>
    sys.exit(main())
  File "/usr/lib/python2.7/site-packages/mapproxy/script/util.py", line 387, in main
    commands[command]['func'](args)
  File "/usr/lib/python2.7/site-packages/mapproxy/script/util.py", line 85, in serve_develop_command
    app = make_wsgi_app(mapproxy_conf, debug=options.debug)
  File "/usr/lib/python2.7/site-packages/mapproxy/wsgiapp.py", line 110, in make_wsgi_app
    services = conf.configured_services()
  File "/usr/lib/python2.7/site-packages/mapproxy/config/loader.py", line 211, in configured_services
    return self.services.services()
  File "/usr/lib/python2.7/site-packages/mapproxy/config/loader.py", line 1770, in services
    new_services = creator(service_conf or {})
  File "/usr/lib/python2.7/site-packages/mapproxy/config/loader.py", line 1857, in wms_service
    tile_layers = self.tile_layers(conf)
  File "/usr/lib/python2.7/site-packages/mapproxy/config/loader.py", line 1789, in tile_layers
    for tile_layer in layer_conf.tile_layers(grid_name_as_path=use_grid_names):
  File "/usr/lib/python2.7/site-packages/mapproxy/util/py.py", line 78, in wrapper
    cache[key] = func(self, *args, **kwargs)
  File "/usr/lib/python2.7/site-packages/mapproxy/config/loader.py", line 1698, in tile_layers
    for grid, extent, cache_source in self.context.caches[cache_name].caches():
  File "/usr/lib/python2.7/site-packages/mapproxy/util/py.py", line 78, in wrapper
    cache[key] = func(self, *args, **kwargs)
  File "/usr/lib/python2.7/site-packages/mapproxy/config/loader.py", line 1445, in caches
    request_format=request_format,
  File "/usr/lib/python2.7/site-packages/mapproxy/config/loader.py", line 1339, in _sources_for_grid
    source = source_conf.source({'format': request_format})
  File "/usr/lib/python2.7/site-packages/mapproxy/config/loader.py", line 747, in source
    coverage = self.coverage()
  File "/usr/lib/python2.7/site-packages/mapproxy/util/py.py", line 78, in wrapper
    cache[key] = func(self, *args, **kwargs)
  File "/usr/lib/python2.7/site-packages/mapproxy/config/loader.py", line 545, in coverage
    return load_coverage(self.conf['coverage'])
  File "/usr/lib/python2.7/site-packages/mapproxy/config/coverage.py", line 107, in load_coverage
    return coverage(geom or bbox, SRS(srs), clip=clip)
  File "/usr/lib/python2.7/site-packages/mapproxy/srs.py", line 89, in SRS
    srs = _SRS(srs_code)
  File "/usr/lib/python2.7/site-packages/mapproxy/srs.py", line 122, in __init__
    self.proj = Proj(init='epsg:%d' % epsg_num)
  File "/usr/lib/python2.7/site-packages/mapproxy/proj.py", line 158, in __init__
    (proj_def, init, libproj.pj_strerrno(errno)))
mapproxy.proj.ProjInitError: error initializing Proj(proj_def=None, init='epsg:28992'): no system list, errno: 2

I didn't have to do any specific configurations in Debian, I simply followed the instructions on MapProxy installation website and everything worked fine but with CentOS it looks like it requires some more settings.

Please can I get help with how this could be solved?


texts: how to extract elevation data using a DEM

texts: I have digitized an LULC layer from satellite imagery, now I want to extract the elevation data from it. I have obtained DEM images from SRTM amd ASTER. How do I proceed?


texts: How do I add a field that determines the feature's transparency in ArcMap?

texts:
I have a feature layer that has points. I want these points to be transparent depending on a field that I've created.

Here are references to being able to actually do this (I do understand that these are mostly tech article and a gis.se question, but somewhere someone has done this):

http://support.esri.com/en/knowledgebase/techarticles/detail/17494
http://support.esri.com/en/knowledgebase/techarticles/detail/21444
Creating symbology for multiple attributes within the same feature class in ArcMap 10

So I believe this is possible, and what I've found is that under the symbology settings, the Advanced dropdown should have a Transparency option. However this is what I find:

Any ideas what I need to enable to allow this?


texts: Unable to export Testing and Training data as TFRecord using Google Earth Engine

texts:
While trying to export training (about 700 points) and testing data (about 300 points) as TFRecord, I am getting the following error:
User memory limit has exceeded

The code I have written, is similar to the material available online, and is:

var band = [ 'B1', 'B2', 'B3','B4','B5','B6','B7', 'B8', 'B9', 'B10', 'B11', 'B12'];
var training = cloudFree.select(band).sampleRegions({
    collection: trainFinal,
    properties: ['CrID'],
    scale: 30,
       tileScale: 16
   });

//Exporting Training and Testing datasets

Export.table.toCloudStorage({
     collection:training,
     description:'Training_Export',
     fileNamePrefix: 'trainingTestData',
     bucket: gcs_bucket,
     fileFormat:'TFRecord'
     });

Is there any mistake in the above code, or should I reduce the number of points / bands?


texts: Plotting points using coordinates from excel

texts: I have an excel spreadsheet with co-ordinates (eg. S31 07 35.3). I would like to use these co-ordinates in arcmap (version 10.2) to plot points based on these coordinates. How do I go about doing this because I cannot use the co-ordinates as is? So basically, I need a way of converting these co-ordinates into something arcmap will understand.


texts:
I am now looking for a code in VBA that can call Tools from ArcToolBox(The little red box)

I have find a Arcobject code(I thought it's can run under VBA) but it's don't work!

I have to delete the URL here....

in the end of the page it's have Arcobject code

But it's don't work!

can anyone tells me how to use it?

thanks for the help!!

@ Mike Toews

Compile error:

Sub or Function not defined

at OpenPointDataset

@artwork21 I call it from Tools>marcos>visual basic editor
It's a simple compiler in GIS not a ui control

I have tried the method with the page http://edndoc.esri.com/arcobjects/9.0/Samples/Geoprocessing/RunToolWithGPDispatch.htm

and I try my Kriging Tool

Sub RunGPTool()
 'Create the Geoprocessor object
   Dim GP As Object
   Set GP = CreateObject("esriGeoprocessing.GpDispatch.1")

  'Set the toolbox
   GP.Toolbox = "C:\Program Files\ArcGIS\ArcToolBox\Toolboxes\Spatial Analyst Tools.tbx"

   'Execute tools
  GP.Kriging "C:\sample.shp", _
    "F30", "C:\AA\456.shp", "250", "12"

 End Sub

but it's highlight error at
         GP.Kriging "C:\sample.shp", _
           "F30", "C:\AA\456.shp", "250", "12"

and says
automation error
unspecified error

how to fix it??

update 2011/06/11 10:00PM(+8)

@Casey

http://edndoc.esri.com/arcobjects/9.0/Samples/Geoprocessing/RunToolWithGPDispatch.htm

in this page

there is no "_analysis"in the sample code...

and I tried GP.Kriging_sa

but it's still don't work....

thanks for you support!

by the way I have tried buffer it's works!!


texts: Error issues when using Cut Fill operation in ArcGis 10.2

texts:
I wonder if someone is able to assist with an issue I have about a series of errors I get when using the cut/fill tool in spatial analyst, ArcGis 10.2? Someone posted a question sometime last year entitled:

https://gis.stackexchange.com/questions/58460/what-is-causing-arcmap-errors-010302-010067-while-using-3d-analyst-cutfill

but as you can see it was closed for being too localized. This, however, is my exact problem, and I have seen that the page has been viewed many times, so I don't think I am the only one! The author of the question said that a miraculous solution to the problem manifested itself the following day; I have not been so lucky. Does someone know of a solution to this?


texts: Clipping vs Joining Polygons to Polygons, which sum of Area should I trust?

texts:
I have two polygon layers, layer 1 and layer 2. I want to know the sum of the areas of layer 2 that are completely inside of layer 1 to include the sum of the intersecting portions of areas from layer 2.  Each item in layer 1 is an area that is larger than each area in layer 2. There are many more items in layer 2.

I first used Clip on layer 2, with layer 1. I calculated the sum of the shape area of the output.

Second I joined layer 1 to the Clip output and selected the sum option.

The summed area of the join output, however, does not match the clip.

Why do these not match? Which should I trust?


texts: After applying the Diff operator of the raster calculator in arcgis to view changes in land use, would like to see the areas of increased and areas of decreased. This is possible using the raster calculator?


texts: Adjusting GPS collected data to a polyline shapefile

texts:
I am collecting data from a GPS that is attached to a car. When I display this data in a map, the set of collected points are supposed to match with the street shapefile. However what I see in my map is points overlying places where the car can't go, Ex. Buildings, Parcels etc.

I understand this is a GPS accuracy issue. My question is: Is there any script or any way to spatially adjust this points to place them on top of the line so the collected data display on its right position in the map.


texts: Why don't Landsat 8 image extents match WRS2 polygons?

texts:
The question doesn't need too many explanations, but the following picture illustrate the problem very well.

I sort of understand why they don't match on width. The satellite drifts away of its path and once in a while it is brought back on the center of the WRS2 path.

On the row though, the image is just cut by the processing software on the ground, so why don't they cut it to match exactly the WRS2 rows?

How do they take the decision of where to cut? For me it would have made sense to cut according to the WRS2, but they are obviously not doing this.

The way I noticed this problem is that the 2 images corresponding to 181028 and 182029 should not overlap, but they actually do a bit.

This image is georeferenced in UTM 35N, which corresponds to the projection of the 2 images.

If you want to recreate exactly the image above, you can download the WRS2 shapefile and the 2 images from the following websites and load them in QGIS.

https://landsat.usgs.gov/pathrow-shapefiles
https://earthexplorer.usgs.gov/
LC08_L1TP_182029_20160805_20170322_01_T1
LC08_L1TP_181028_20160814_20180527_01_T1


texts: How do I stack feature class annotation?

texts:
I have map annotation that is stored in a database (as a feature class), but that is not linked to the features. For the life of me, I can't figure out how to edit it so that I can take a label of say 3 words and stack it so that one word is on top of the other, like so:

From:
How To Stack

To:
How
To
Stack

Any help would be greatly appreciated!


texts: Exporting list of values into csv or txt file using ArcPy?

texts:
I would like to know how to export in ArcGIS Desktop a list of values calculated in Python script into one of the following data formats: csv, txt, xls, dbase or other. I would also like to know how to create such file in case that it doesnt exist.

The list of values looks like res=(1,2,3,...,x). Each value must be written into a new row.


texts: Strahler calculation ignoring large rivers

texts: I have a DEM of a portion of the Yukon. I merged numerous DEM rasters then selected a subset for my AOI. I filled the sinks and then calculated the Strahler, but it seems to be ignoring large streams and rivers. The image shows Strahler >4 for reference.


texts: Creating a Service Area with 6 hours Travel Time on ArcGIS Online

texts:
I am creating a series of Service Areas using the built in networking features on ArcGIS Online. I need to create an area with 1 hour intervals for 6 hours of travel time. The issue is that ArcGIS Online does not allow for any longer than 5 hours through the online networking tool, and unfortunately I need that 6th hour!

I realize I may be out of options for accurate analysis, so I have moved on to other methods to try and at least show an approximate 6th hour. Buffering the largest area does not work as it returns a smooth area that just doesn't look good at all. I have tried creating a copy of the 5th hour, and using the Stretch Proportionally tool in my edit vertices, but this does not work well either.

I was hoping to get a couple more ideas I could try, because I am not sure what other options I have.


texts: Export color corrected rasters from a catalog in arcmap

texts:
I have several raster in a raster catalog for which I(ve used color correction. I'de like to export each raster as a geotiff with exactly the same resolution and world file than the original geotiff rasters. By exporting the selected rasters from the raster catalog with colors corrected, in a blank raster catalog (MANAGED), the image file is not created.
Can you explain how to export color corrected rasters to georeferenced image files with teh same properties than the input rasters ?


texts: How to animated feature marker in openlayers like in google map example of bouncing marker. I'm new in openlayers so if you have a sample code much better.


texts: Is a new map object required in order to create a wms overlay on top of a current map object with wfs layers using the leaflet api

texts: I am trying to add a wms tile layer overlay over a wfs tile layer which I added to a map object. When I tried to simply add a layer on top of my current setView for my map object it overwrote the first layer's and all the GUI I added was gone. How do I add a wms tile layer overlay?


texts: Adding information to point when created with PyQGIS

texts:
How do I add information to a point when I'm creating it with Python in QGIS. I want to add some info to the point like id, name, etc...

This is my code:

layer = QgsVectorLayer('Point?crs=epsg:4326', 'point' , 'memory')
prov = layer.dataProvider()
point = QgsPoint(point)
feat = QgsFeature()
feat.setGeometry(QgsGeometry.fromPoint(point))
prov.addFeatures([feat])
layer.updateExtents()
QgsMapLayerRegistry.instance().addMapLayers([layer])


texts: Change Layers Panel background color in QGIS

texts: In QGIS it's easy to change background color of the map canvas, but is there a way to do the same for the Layers Panel ? In cases of very light colored styles these appear not so clear (if at all) on a white background. I know that in QGIS 3 there is an option to toggle GUI theme from default to Night Mapping, wich turns anything into darkness, but is there a way to set Layers Panel background to e.g. black independantly from the GUI theme? Python perhaps?


texts: Trigger to fill attributes (data, client ip) in table for each new objects, PostGIS

texts:
In PostgreSQL/PostGIS database stores shapefiles. Client users software - QGIS. Please help me write PostgreSQL/PostGIS trigger, which for each newly added/updated feature in layer added 'time create' and 'client ip address' in attribute table.
Attribute table contain next fields: 'gid' - primary key, 'updated' - data create/update feature, 'nameuser' - client ip address.
I write next trigger:

CREATE OR REPLACE FUNCTION process_audit() RETURNS TRIGGER AS $audit_roads$
BEGIN

IF (TG_OP = 'INSERT') THEN
 INSERT INTO roads SELECT now(), inet_client_addr(), NEW.*;
 RETURN NEW;

 ELSIF (TG_OP = 'UPDATE') THEN
 INSERT INTO roads SELECT now(), inet_client_addr(), NEW.*;
 RETURN NEW;

 END IF;

 RETURN NULL;
 END;
$audit_roads$ LANGUAGE plpgsql;

CREATE TRIGGER audit_roads
AFTER INSERT OR UPDATE OR DELETE ON roads FOR EACH ROW EXECUTE PROCEDURE process_audit();

but when I add new feature in QGIS and save them, I receive message: "INSERT has more expressions than target columns"

Tracking Edit History using Triggers - there is example of Tracking Edit History using Triggers, but history saved in nyc_streets_history table (not in nyc_streets)

Thanks!


texts:
I have a number of vector layers, I want to be able to drag a feature off the map and drop it into a native application.  Note, I do not want to move the features around, I already do that when in a particular mode.

I've already got the transfer working from a list I build up in the page when users select features on the map using html5 drag and drop, but it's a bit clunky having to select the feature then drag and drop the list item.

I can not see anything in Openlayers that will allow me to handle the ondragstart event for a feature :-(


texts: Exporting ArcGIS TIN file to *.tin file?

texts:
How can I create a .tin file from a ArcGIS TIN file?

I can't seem to find a export to .tin solution.


texts: Web AppBuilder Widget Configuration Validation

texts:
I have created a custom widget with configuration in web appbuilder. What I need to do is disable "Ok" button or do nothing if all fields in the config screen are not filled. In other words, I need to do save only when all inputs are filled. Any idea how to achieve it?

Regards,
Jamil


texts:
I'm new to GIS and mapping, and I need to learn how to create web map application for my project. I'm really confused on what to use. I have read PostGIS, OpenStreetMap, and OpenLayers. I'm wondering how those API's work together. I'm an experienced PHP programmer, so is there a way that I can apply my PHP skills?

Later on this project, I'll be calculating the shortest path and creating a routing guide. Can you guys suggest a best guide for accomplishing my project?

Also, is there a way that I can use only domain server or host like bluehost? Since I'm only required to work within a small town.


texts: How to add BC-dates (age of a continent)

texts:
I want to simulate continental drifts using the time tool of QGIS (V. 3.16.0). Turns out I can only insert certain dates into my attribute table. I can't even type in the date 0000-01-01 without getting a NULL-field. Typing in a negative value is just not possible.
I cannot get QGIS to read a BC or AD after my date, as suggested in a response to this question.
It would be okay for me to define the year of when the landmasses formed to the year 0000. But I cannot define years with more than 4 digits without simply repeating the date! (yyyy → 2020; yyyyyyyy → 20202020 instead of 00002020)
So what method would you suggest to define how many million years before Christ a feature began to exist in the current form?


texts: Image.constant: Parameter 'value' is required. Error in Google Earth Engine

texts:
I'm trying to create the Vegetation Condition Index with Landsat time series images in Earth Engine, but I get the error:

  ImageCollection (Error) Error in map(ID=LC08_198032_20130813):
  Image.constant: Parameter 'value' is required.

Below is the script code, maybe something I was doing wrong at the section "//Create VCI index":

var l8collection = l8.filterDate('2013-04-01','2019-11-30').filterBounds(roi);
print(l8collection);

// Mask pixels with clouds and cloud shadows using the 'pixel_qa' band

var maskClouds = function(image){
// make a new single band image from the pixel qa band
  var pixel_qa = image.select('pixel_qa');
  // retain clear (0) and water (1) pixels
  return image.updateMask(pixel_qa.eq(322));
};

// use "map" to apply the function to each image in the collection
var l8masked = l8collection.map(maskClouds);
print(l8masked);

// create function to add NDVI using NIR (B5) and the red band (B4)
var getNDVI = function(img1){
  return img1.addBands(img1.normalizedDifference(['B5','B4']).rename('NDVI'));
};

// map over image collection
var l8ndvi = l8masked.map(getNDVI);

// filter images not containing data in the study area
var filtered = l8ndvi.filter(ee.Filter.neq('count', 0));

print(filtered);

// Create VCI index
var vci = filtered.map(function(img){
 var id = img.id();
  var min =  img.reduceRegion({reducer: ee.Reducer.min(),geometry:roi, bestEffort:true, scale: 30}).get('NDVI');
 var max =  img.reduceRegion({reducer: ee.Reducer.max(),geometry:roi, bestEffort:true, scale: 30}).get('NDVI');
 var Index = img.expression("(NDVI-min)/(max-min)",
 {"NDVI" : img,"max" : ee.Number(max),"min" : ee.Number(min),})
 .rename('TCI')
 .copyProperties(img,['system:time_start','system:time_end']);
 return Index;
});
print(vci);


texts: Implementing custom (Vincenty/Haversine formula) functions for MySQL vs Using PostGIS?

texts:
I currently am using MySQL for all of my relational data. I have hit a crossroads with geo calculations. Since, MySQL does not have good support for spatial functions, I am contemplating on using Postgres's PostGIS extention.

But, for this I will have top run 2 DB servers (MySQL and Postgres) as I do not want to move away from InnoDb(for all the goods that it offers).

Does implementing the functions cause any performance issues(lag)?
Is the trouble worth it or should I move ahead with PostGIS?
Is moving the whole DB from MySQL to Postgres any good (Although I am not exactly considering this)?

Also, is implementing Haversine or Vincenty better? I do understand that Vincenty formula provides much better accuracy. But is there too much cost on CPU?


texts:
I am using the following code to plot weather readings on a map.  It works but there is a very noticeable and annoying blink each time new readings are plotted.

function draw_map(current_obs) {
  vectorSource.clear();
  for (stn in current_obs) {
    addMarker(current_obs[stn]); // does vectorsource.addfeature()
  }
  vectorSource.changed(); // needed to redraw after adding markers
}
i=0;
//time_series is an array of current_obs data
loop_timer = setInterval(() => {
  draw_map(time_series[i]);
  i = i + 1;
  if (i >= timeseries.length) { i = 0 };
}, 500);

Is there a way to make the transition between time steps visually smoother? I've thought about doing features individually using removeFeature()/addFeature()/changed() but I don't how to map existing features to station observations.  (There are about 40 stations in the dataset, each having several parameters plotted on a single feature)

I'm using OpenLayers 5.3.0

Following the suggestion below I tried re-writing draw_map()

function draw_map(current_obs) {
features = []
for (let stn in current_obs) {
    //make_Marker returns a feature object rather than adding it
    features.push(make_Marker(current_obs[stn]))

}
vectorSource.clear();
vectorSource.addFeatures(features);
vectorSource.changed(); // needed to redraw after adding markers

}

However this did not seem to make a noticeable difference.  So I tried refreshing each station individually (easy once I realized you can  add arbitrary properties to a feature object).

function draw_map(current_obs) {
for (let stn in current_obs) {
    new_feature = make_Marker(current_obs[stn])
    vectorSource.forEachFeature(
            function(ft) {
                foundit = (stn === ft.get('station_name'))
                if (foundit) { vectorSource.removeFeature(ft) }
                return foundit
            })
    vectorSource.addFeature(new_feature);
    //new_feature.changed(); //this didnt seem to be needed
}

However even this did not seem to reduce the blinking.

My markers are using an SVG path for the ol.style.Icon element (a wind direction arrow) and commenting that out seems to greatly reduce the blinking so perhaps its just that rendering SVG takes too long.  I will have to re-think that part.


texts: Creating raster from grid algebra expression using point attribute values?

texts:
I have a point shapefile representing groundwater production wells with annual pumpage values distributed over a general area. I'd like to create a raster layer with each cell representing the potential production impact weighted by distance using the following expression:

outputraster = SUM n/i (Xi/log10(Di))

Where 'Xi' is the annual production of well i and 'Di' is the distance in meters from well i.

Since there a multiple wells scattered over a wide area, the resulting raster cell values should be highest based on the density of wells in a given area, proximity to a given well, as well as the production value of a given well.

I've begun by creating a raster dataset, but am having trouble initiating the map algebra expression.

Can this calculation be performed within the raster calculator tool?

I have very little coding experience and my work with rasters extends only to interpolation and basic raster algebra expressions.


texts: Comparing multiple field values of two feature classes using ArcPy?

texts:
I incorporated the following code into my script from Comparing field values of two feature classes using ArcPy? and it works great but I also need to compare two other fields values (i.e. YEAR & MONTH) from the same two FC's.

There is no problem with the code the way it is, I just need help on modifying the code to check 3 fields' attributes from 2 FC's to see if each row for the three fields contain the same attributes and if so not to copy the row to the target feature class. The fields in each FC do match exactly.

import arcpy

#input feature class
fc2 = r"C:\temp\temp.gdb\fc2"
#target feature class
fc3 = r"C:\temp\temp.gdb\fc3"

#check field
checkField = "SRNumber"

#Get list of values in field from target feature class
checkValues = [r[0] for r in arcpy.da.SearchCursor (fc3, checkField)]

#Get list of fields
fields = [f.name for f in arcpy.ListFields (fc2)]

#Get index of check field
index = fields.index (checkField)

#Create insert cursor for fc3 to allow appending of rows
inCursor = arcpy.da.InsertCursor (fc3, fields)

#Create search cursor to iterate input feature class
cursor = arcpy.da.SearchCursor (fc2, fields)
#iterate
for row in cursor:
    #Get value to check
    checkValue = row[index]

    #Check if value is in target fc
    if checkValue in checkValues:
        #skip if value is in target fc
        continue

    #Insert row otherwise
    inCursor.insertRow(row)

#Clean up
del cursor
del inCursor


texts:
I am trying to create a transit map with certain details such as landmarks and important centers for reference in it.

I have mapped the lines and stops ion QGIS (2.16.2). The base map is Google street maps.

What I want to do is take this whole project as a vector image, and place it somehow in Illustrator, so I can work with the map as a reference to design the map.

What this will allow me to do is use the simplest property of a vector image- the ability to zoom in to parts and view in detail. Not only will I be zooming in the AI document, but also the map so I can view and edit in parallel.

Currently, I am using  a raster image which pixelates a lot when I zoom in.

Any idea how?


texts: How to calculate the area within a 4m clipped buffer derived from a raster

texts: I am using QGIS 2.8.3. I started by classifying my area into habitat types using the semi-automatic classification plugin (SCP). Then I added points, and a 4m buffer around the points. So I could then clip the results of the SCP by the 4m buffer, to figure out what the habitat types are within the 4m buffer. Now I want to know the area or percentage cover of each habitat type for each point within the 4m buffer based on the SCP classification. Any ideas on the best way to do this?


texts: How to efficiently (programmatically) determine whether two datasets have / don't have identical content?

texts:
Given:

a geodatabase A that contains 100 datasets named a1 .. a100
a geodatabase B that contains 100 datasets named b1 .. b100

I want to programmatically (*) determine for each dataset pair (ai, bi) whether they have identical content. Since I am comparing 100 pairs, I need an efficient comparison method. Ideally, the whole comparison would take only a few seconds.

(*) Note: I mention the term "programmatical", not because I am looking for code examples (though I'd gladly accept them), but in order to emphasize that I am looking for a very speedy comparison method, which would never be possible when comparing 100 dataset pairs manually.
I am planning to implement this comparison method myself, so what I am essentially looking for is an algorithm, and not a ready-to-use tool (unless perhaps when it's open-source).

I am aware that I am likely asking the impossible, since this would require comparing the datasets' complete contents (perhaps with the tools in the Data Management → Data Comparison toolset); or at least comparing dataset hashes / digests, but generating digests would also require going through all of the datasets' data first.
Therefore my best approach so far is the following:

Determine first which dataset pairs (ai, bi) cannot possibly have identical content.

Perform a full data comparison only for the remaining dataset pairs.

My questions:

Does ArcGIS happen to auto-compute some kind of dataset digest that I could query? If so, how?
(I am not aware of anything of that sort, so I expect the answer to be "no". Please prove me wrong.)

What are some very efficient, reliable ways of determining whether two datasets cannot possibly have identical content?
(I have so far considered comparing modification timestamps, though I don't know how reliable these are, and comparing the datasets' schemas. How reliable are timestamps in an ArcGIS geodatabase? Are there other dataset characteristics that could serve for this purpose?)


texts: Showing just some labels based on feature intersection of layers using QGIS

texts:
I'm trying to use "expression based" labeling in QGIS, but I'm stuck with a particular case.

I managed to use it for a project with the help of other topics and FAQ but I'm not skilled enough to take care of this problem alone.

I have 5 layers, one "project" layer with around 10 objects (10 different marinas around rivers in France) and 4 layers for particular water area (different scales) with also polygons. I want to only show labels of those 4 last layers of the polygons that intersect any of my 10 objects of layer 1.

What query do I need to use?

Note: I already used the vector selection option and created a new layer with the selection. I'm looking for another option, label based because this is too long. I know there is also a "show/hide" option, but it is not very simple to use in my opinion.


texts: R 'gstat - Warning: singular model in variogram fit

texts:
I am trying to perform a kriging interpolation. But when it comes to the point when i want to fit the variogram to the model I get the following warning:
In fit.variogram(vario, vgm("Sph")) : singular model in variogram fit

I am not sure where I have to look at to solve the porblem.

my data (these are just 10 / 216 lines):

Place;Latitude;Longitude;Temperature;Humidity;Windspeed;AirPressure
Aachen;50.77999878;6.09999990;3;93;15;1020
Abbikenhausen;53.52999878;8.00000000;7.9;83;1.9;1022
Adelbach;49.04000092;9.76000023;3.1;91;8.0;1014
Adendorf;51.61999893;11.69999981;1.9;76;2.9;1018
Alberzell;48.45999908;11.34000015;4.6;97;1.9;1012
Altenstadt;47.83000183;10.86999989;3.8;90;3.1;1012
Altersteeg;51.58000183;6.32000017;4.4;89;1.0;1017.5
Angermuende;53.02999878;14.00000000;1.5;89;1.0;1020
Arnsberg;51.11999893;7.32999992;2.3;100;13.0;1018

my code:

library(raster)
library(sp)
library(gstat)

WU_data <- read.csv(file = "./WU_Data.csv", header = TRUE, sep = ";")
WU_data <- WU_data[complete.cases(WU_data),]

min_lon <- min(WU_data$Longitude)
max_lon <- max(WU_data$Longitude)
min_lat <- min(WU_data$Latitude)
max_lat <- max(WU_data$Latitude)

Longitude.range <- as.numeric(c(min_lon,max_lon))
Latitude.range <- as.numeric(c(min_lat,max_lat))

grd <- expand.grid(Longitude = seq(from = Longitude.range[1], to = Longitude.range[2], by = 0.1),
                   Latitude = seq(from = Latitude.range[1],to = Latitude.range[2], by = 0.1))  # expand points to grid
coordinates(grd) <- ~Longitude + Latitude
gridded(grd) <- TRUE

WU_data_spatial <- WU_data
coordinates(WU_data_spatial) = ~Longitude + Latitude

vario <-  variogram(Temperature ~1, WU_data_spatial)
vario.fit <- fit.variogram(vario, vgm("Sph"))

plot(vario):


texts: Can QGIS automatically load an SLD with a Shapefile?

texts:
I have SLD's/lyr files with my shp/gdb files. They all have the same name except for the extension. Is it possible for me to force qgis to load the sld with the shp so that I don't need to manually go to properties, load style and ok it?

In arcmap I can just open the lyr instead of the .shp -is there a similar process for qgis?


texts: Can I have the "store relative pathnames to data sources" box checked as the default?

texts: Is there a way to have the "store relative pathnames to data sources" box checked as the default rather than always having to go into the document properties and check the box? ArcMap 10.0.


texts: How can two Oracle Spatial shapes be disjointed and interacted at the same time?

texts:
This Oracle Spatial query indicates a polygon and a line have interaction but also that they are disjointed.
How can that be?

--disjointed but interact???
SELECT lines.pavement_key, poly.SUB_DISTRICT_NAME_UPPER ,
   SDO_GEOM.RELATE(lines.shape, 'DETERMINE'  , poly.shape,.5) detemination,
   SDO_GEOM.RELATE(lines.shape, 'ANYINTERACT', poly.shape,.5) ANYINTERACT
FROM DW_RAH_PAVEMENT_SECTIONS_DIM lines join DW_SUBDISTRICT_MASTER_DIM poly
on MDSYS.SDO_ANYINTERACT(lines.shape, poly.shape)='TRUE'
where POLY.SUB_DISTRICT_NAME_UPPER='CRAWFORDSVILLE'
and LINES.PAVEMENT_KEY in (1003);

below is a map image which shows the line is actually in the polygon


texts:
I am trying to use python to parse the .pvl metadata file from some GeoEye imagery. I would like to write a script that locates the productCloudCover and countryCode tags found within the .pvl. Does anyone know where I could find information on how to do this?

The structure of the file is as follows.

copyright = "######"
companyAddress = "########"
BEGIN_GROUP = "product order info"
     projectName = "#######"
     BEGIN_GROUP = "#####"
         productCloudCover = "0.2"
     END_GROUP
END_GROUP = "product order inf"

Thanks