texts: How to profile rendering performance using QGIS 3.10.x or above?

texts:
I am looking to be able to profile the rendering performance of a map and/or its individual component layers using QGIS.

The metrics I am most interested in are:
Access time
Storage IOPS
Storage time
Network IOPS
Network time
Rendering time
Rendering perf (features/sec or tiles/pixels/sec)

What I am attempting to determine are best-practices for simplification and storage of my data to ensure that it can be viewed and interacted with in the most responsive manner.

A few for-instance questions I'm looking to gain insight on using profiling:
Does FGBD serve faster than SHP?
GeoPackage vs FGDB/SHP?
Does GPKG/SQLite3 page size make a difference?
Does subdivision to fit given geometry segments in a single database page make a difference?
Should I index X/Y/Z attributes in the table, or does the increased database size mitigate the speed gain?
Is it faster to have a geometry table joined to an attribute table, or one aggregated table of geom & attributes both?

I've been researching numerous methods of optimization, but it is hard to find concrete methods of profiling these against one another on my target machine running QGIS.

I currently have a Processing Model that does a few "optimization" steps for me against any input data I have, but I would certainly be open to expanding its repertoire and standardizing upon these best practices.

Currently, I have Snail installed, but it doesn't show nearly as much as the Task Manager that ships with W10 does, so it hasn't been of much use to me for profiling.

Similarly, messages in the Log Messages panel aren't as helpful as I would like.


texts: Getting a count of points in a multipolygon?

texts:
I have a set of lottery retail locations, latitude and longitude, that I transformed to points and added the geometry column like this:
ALTER TABLE lottery.retailers ADD COLUMN geom geometry(Point, 4326)
UPDATE lottery.retailers SET geom = ST_SetSRID(ST_MakePoint(retailers.lon, retailers.lat), 4326)

That seemed to work well.
I also have a set of polygons that are the county boundaries of Iowa, USA. Both tables are set to SRID 4326.
https://geodata.iowa.gov/dataset/county-boundaries-iowa
I am trying to count the number of retail locations in each county with this query, but it just gives me a count of 0 for each county listed:
Select county.county, count(retailers.retailer_id) as total
from lottery.county
left join lottery.retailers on st_contains(county.geom, retailers.geom)
group by county.county

Obviously, I am doing something wrong here, but it doesn't appear to be far off. I feared that the coordinate points aren't right because I transformed them from addresses, but I double-checked 10-15 of them and they were all right. Any idea on how I can adjust this query to get the count?


texts:
I have a GeoTIFF which is at 30m resolution. I need to downsample it to 250m resolution (using some aggregating function such as mean, mode, sum). My stack is python (rasterio and GDAL bindings)

I can find a way to do this directly using numpy writing some nasty double loops and iterating over the matrices directly but looking for a cleaner way.

Note that the main issue is that 250 is not a multiple of 30 and there will be some overflow. For instance if the aggregating function is the sum the algorithm would be to use:

8 X 8 block which is 240m X 240m in base resolution
Chose a partial pixel strip at 10m X 10m to complete the sum to 250m
The value for 2. would be area weighted and added to the sum from 1.

Please let me know if I was unclear. Again I can find a brute force way to do it using matrix manipulations but was wondering if there is an elegant way. I am open to using R also and read about the projectRaster and aggregate methods but they seem to work on integer cut offs also.

In general it may be really cumbersome to do this for an arbitrary resolution say 31.567m or so.


texts:
I am trying to create a map of Presidential election results by ZIP code, but all I have been able to find is election data by county (too broad) and precinct (too narrow).
Does anyone know of a public data set that maps election data to ZIP code?


texts:
I have question -i need very easy shift 1 point. It is very difficult in new version.
1.Why I can't shift point closer like red distance? Between red and yellow is not possible snapp. Blue is what I need. I had try change snapping options tolerance, radius...but same result. Version is 3.12.

why I cant shift point with left mouse continually? I need click and release buttton -shift and then click again.Terrible. In old version I could easy chatch point and move, here not.
I upload short video here, you can look:
https://uloz.to/file/N2X3kSbKYzi0/20200408-104927-mp4

I dont know why yesterday closed my question regarding this?It is very serious problem for all I think. If not solve i need go back to old Qgis. Because It is very comfortable I need catch and shift every point excatly where i need every milimeter.
Thanks]1


texts: Error code 500, error performing execute operation from Print task of ArcGIS API for JavaScript?

texts:
I'm trying to use print task to generate a PNG or PDF file of my web map but i always get an Error code 500, error performing execute operation.

Any Ideas?

I'm using ArcGIS API for JavaScript 3.2.0

var dojoConfig = { async : true };

on(dom.byId("png"),'click', function () {
printTaskSetup("MAP_ONLY", "png32");
})
on(dom.byId("landscape"),'click', function () {
printTaskSetup("A4 Landscape", "PDF");
})
on(dom.byId("portrait"),'click', function () {
printTaskSetup("A4 Portrait", "PDF");
})

function printTaskSetup(layout, format){
var template = new esri.tasks.PrintTemplate(
template.format = format;
template.layout = layout;
template.layoutOptions = {
 "scalebarUnit": "Miles",
 "copyrightText": "Clerk",
 "showAttribution": false
}
template.preserveScale = true;
var params = new esri.tasks.PrintParameters();
params.map = map;
params.template = template;
printTask = new esri.tasks.PrintTask(print URL, dojoConfig);
printTask.submitJob(params, printResult, printError);
}
function printResult(result) {
  window.open(result.url)
}
function printError(error) {
  alert(error);
}


texts: Is there a way to rename relationship class labels (forward and Backward labels)  after I've created them ? or I have to delete them and create them again ?


texts: Proj4 projection supports WebGLPointsLayer but throws GeoJSON polygons off

texts:
In order for the following WebGLPointsLayer

    this.layer = new WebGLPointsLayer({
      source: new VectorSource({
        url: "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_month.geojson",
        format: new GeoJSON(),
        attributions: 'USGS'
      }),
   })

to render within the correct extent in a map that has a view center shifted to [103.8198, 0], I am configuring the following custom proj4 projection

proj4.defs('sphmerc180', '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=103.8198 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs');
register(proj4);

var sphmerc180 = getProjection('sphmerc180');
sphmerc180.setExtent(getProjection('EPSG:3857').getExtent());
sphmerc180.setGlobal(true);

and then writing the view like so:

this.map = new Map({
  target: 'map',
  layers: this.layers,
  view: new View({
    projection: 'sphmerc180',
    center: fromLonLat([103.8198, 0], 'sphmerc180'),
    zoom: 1
  })
});

This works nicely to shift the WebGLPointsLayer from this extent:

into the expected extent:

However, this technique wreaks havoc on this GeoJSON polygon layer I'm using from Natural Earth:

which otherwise renders normally with the center shift to [103.8198, 0] without the custom proj4 projection. I understand that some polygon coordinates are being sent round the world, which is why many polygons are rendering that way.

What is the best way to accommodate both of these layers?


texts: Problems closing dbf-table edited with ArcObjects

texts:
To change values and add new values to a dbf-file, I connect to the dbf and use the ITable interface to make the edits. At the end I release the COM-objects of the table and the workspace. But when I try to open the changed file, it says the file is read-only. Did I forget to release something else? Or is there another reason for the read-only status?

Thanks and regards
Fab


texts: Count points in polygon (weighted) in QGIS 3

texts:
In QGIS 2 you could perform a weighted point within a polygon layer using Processing toolbox > Vector Analysis > Count points in polygon (weighted).

I could well be missing something, but I can't find the same function in QGIS 3 (screenshot below). Has it been superseded by another function?

Maybe it requires a plugin?


texts: Is there any documentation / definition of the operations within, crosses, intersects etc that can be found on the spatial query just like ArcGIS provide?


texts:
I have a special map. Its rather a schema not a common world map. Not very high but it is very long. Now when I display the map, I see a part of the schema at zoom level 0 (I additionally set it programatically usin setZoom(0)) also the zoom slider is all the way down and with -(minus) button I cannot zoom out more to see the whole schema.

I set the map like this:

var TSKmap = new MyMap({
    target: 'map',
    view: new ol.View({
      zoom: 1,
      center: [0,0],
      projection: ol.proj.get('EPSG:4326')
    })
});

Add WMS layer using TileWMS source with the following parameters

BBOX    =-90,-180,270,180
CRS     = EPSG:4326
FORMAT  = image/png
HEIGHT  = 256
LAYERS  = lyr1,lyr2,lyr1_labels,common_labels
REQUEST = GetMap
SERVICE = WMS
TRANSPARENT = true
VERSION = 1.3.0
WIDTH = 256
idk   = 17
idp   = 101
map   = schemas/schema1.map
version = 20.3

The l1 layer is setup in the schemas/schema1.map like this:

LAYER
NAME "..."
METADATA
  "wms_title"          "..."
  "wms_abstract"       "..."
  "wms_keywordlist"    "..."
  "wms_extent"         "-5000 -2000 50000 4000"
  "wms_enable_request" "*"
  "wms_include_items"  "all"

  "wfs_typename"       "..."
  "wfs_version"        "1.1.0"
  "wfs_title"          "..."
  "wfs_abstract"       "..."
  "wfs_keywordlist"    "..."
  "gml_include_items"  "all"
  "gml_featureid"      "featid"
  "wfs_enable_request" "*"
END

Interesting information from mapfile:

MAP
  NAME "TSK_WMS_WFS"
  DEBUG 5
  SIZE 1020 690
  EXTENT -5000 -2000 50000 4000
  UNITS METERS
  FONTSET etc/fonts.txt
  MAXSIZE 4096

  WEB
    IMAGEPATH "tmp/"
    IMAGEURL "tmp/"
    METADATA
      ...
      "wms_srs"                           "EPSG:4326"
      "wms_encoding"                      "UTF-8"
      "wms_bbox_extended"                 "true"
      "wms_extent"                        "-5000 -2000 50000 4000"
      "wms_server_version"                "1.3.0"
      "wms_enable_request"                "*"

      "wms_feature_info_mime_type"        "text/html"
      "wms_format"                        "image/png"

      "wfs_srs"                           "EPSG:4326"
      "wfs_enable_request"                "*"
      "wfs_encoding"                      "UTF-8"
      "wfs_extent"                        "-5000 -2000 50000 4000"
      "wfs_version"                       "1.1.0"

    END
  END

  PROJECTION
    "init=epsg:4326"
  END

  ...

Any Ideas how to set the view so, that I shoul see the whole schema when I zoom to 0?


texts: how can i use openlayers to use my locally rendered tiles

texts: I've been stuck on how to use my local tiles which i rendered with mapnik's generate_tiles.py and all the data is on postgresql. I'm kinda confused on what i should do to "feed" the tiles to openlayers. The map would be used in django btw. hoping you could help me with this.


texts: Cell value changes during raster clipping using QGIS?

texts:
I have created a raster with a range of values from -1 to 1. I would like to extract a section by clipping to a polygon extent. I've been using Raster -> Extraction -> Clipper and QGIS has been giving me outputs clipped to the correct extent.

BUT - The data in each cell of the clipped versions of the raster are different from the original. For example, when there should be no data I'm getting numbers like 0.02 and 0.0314. Also, the field for nodata value = 3.40282e+38.

If I use gdalwarp I can change the nodata value to something outside of my -1 to 1 range (like -10 for example), but the cell values in the clipped raster are still different from my original raster. It seems that some kind of transformation is happening to the data during clipping.

Has anyone encountered this before?


texts: QGIS is not showing/rendering points that are present in MySQL server

texts:
I have a MySQL database in which there are auto generated points almost 1000 created via Laravel (PHP framework) it has a very simple schema as shown below. I just want to visualise them on an OpenStreetMap map in QGIS,

I did add a vector layer and set the source to be my points table
The points do not show on the map, but when I open attributes table and then keep moving to next feature from the attributes table window the QGIS blinks red circles on the display indicating that there are points there. But unfortunately it does not render and I have been struggling for a while to get them displayed. I have attached another screen shot showing when it blinks them. I am new to QGIS.

UPDATE :
After doing a lot of search and experiments with changing data inserting it manually in mysql using different variations, turns out this is an OLD QGIS Bug. There is a work around to to make it render the data but then it won't allow you to edit the data.

Here is the link, issue already reported but not yet fixed as I tested in QGIS 3.10 version

https://issues.qgis.org/issues/20563
https://github.com/qgis/QGIS/issues/28383


texts: Raster appears unaligned with other data layers?

texts:
We use the Dallas-Ft Worth rasters for our internal processes. We use publicly-accessible data to superimpose markings representing runways and runway ends (the dark gray line and green dots on the uploaded image, respectively). The coordinates for runway ends, for instance, perfectly match the data found here. But for some reason, these markings don't match up with the FAA map. Not even an .osm file from OpenStreetMap (rendered as the transparent purple airport) matches up with the FAA map.

Why is that?

Is this an issue with how ArcMap renders the Dallas-Ft Worth raster from the FAA?

Is there some sort of configuration we need to tweak in importing the raster into ArcGIS?


texts: How to batch produce a state map with individual counties highlighted?

texts:
I'm working on a basic choropleth map that visualizes one variable for all the counties in a single state. I want to create a map for each county that highlights (or adjusts the border style for) the relevant county boundary. The map will otherwise be identical, showing the whole state, the same variable classes, and so forth. If easier, other style tweaks that draw attention to a specific county would also be suitable. Maps will be used in printed documents.

I'm primarily using QGIS 2.1, but solutions that involve ArcGIS would be fine. Python or other frameworks are acceptable as well, but please reference some documentation if you can so that I'm able to brush up on relevant details.

I've looked into trying to use the QGIS Print Composer Atlas, but it doesn't seem that Atlas alone will yield a solution.

In other words, what's the best way to batch produce maps (to be exported as images) that stylize each individual polygon feature while maintaining the full map extent and other details?


texts: Where to find shp2pgsql windows binaries?

texts:
Is there anywhere available shp2pgsql binaries for windows without installing the whole PostGIS? I particularity need the command line version of that import tool.

I would like to share a batch script that uses the shp2pgsql.exe to my co-workers but don't know how to share the shp2pgsql.exe. Could I just copy that file from my postgres/bin folder? I suppose that then I need at least libpg.dll. How about other files (libgdal-1.dll, libproj.dll)?


texts: How to I recreate MultiPolygon with Rings/Holes using org.geojson

texts:
I have a few original MultiPolygons that contain holes...I would like to extract the Polygons and recreate them WOTHOUT the holes.  how is this done?
Original MultiPolygon.
https://gist.github.com/boundaries-io/978eaa4a10df9467638a5eb9259c84e6

when breaking up the Multipolygon
   org.geojson.MultiPolygon reCreateMultiPolygon = new org.geojson.MultiPolygon();
             for (List<List<LngLatAlt>> polygonPoints : multiPolygon.getCoordinates()) {
                for (List<LngLatAlt> point : polygonPoints) {
                    Polygon p = new Polygon(point);
                    reCreateMultiPolygon.add(p);
                }
            }

             String fileName2 =   saveToFileSystem(reCreateMultiPolygon);

The this is treating each polygon as its own polygon...I am unclear how to recreate this with the holes intact...not treat them as normal polygons.
results:/
Incorrect recreated MultiPolygon https://gist.github.com/boundaries-io/89e51458e9f5509b5457f9033474bc2d

my end goal is to create a convert that does a org.geojson.MultiPolygon to rg.springframework.data.mongodb.core.geo.GeoJsonMultiPolygon   to save in Mongodb via spring data mongodb.   if this can be done via JTS,geotools,etc. thats fine...


texts: Choosing weight for red, blue and green using Pan-Sharpened-Raster-Dataset of ArcMap

texts:
I am pan-sharpening satelite image using the ArcMap tool Pan-Sharpened-Raster-Dataset.
What weight should I use for red, blue and green and how I can decide about that?


texts:
I am trying to use the SEXTANTE module GRASS (and SAGA) for vector analysis but whatever function I use I always get the same error:"Could not load layer: G:/closest.shp
Check the SEXTANTE log to look for errors in algorithm execution"

SEXTANTE log says (for "v.distance" module):
"Sextante.runalg("grass:v.distance","G:/For QGIS/example_1/point.shp","G:/For QGIS/example_1/point_II.shp",1,"col_name","id",False,False,"G:/closest.shp")"

I get this on three different computers with QGIS 1.8.0., SEXTANTE 1.0.7., updated Win 7.

I think this question is the similar problem:
https://gis.stackexchange.com/questions/24895/qgis-sextante-doesnt-create-load-output


texts: How to create feature LineString in a Polygon shapefile layer (OGR/GDAL)?

texts:
I used below codes to create shapefile layer:

OGRDataSource       *poDS,  *poSHPDS;
OGRLayer   *poSHPLayer;
...
poSHPLayer = poSHPDS->CreateLayer (shpFileName, NULL, wkbPolygon , NULL);

And I used below code to create feature:

OGRFeature *poSHPFeature;
poSHPFeature = OGRFeature::CreateFeature (poSHPLayer->GetLayerDefn());
poSHPFeature->SetGeometry (poGeometry);

But when I tried to add create a LineString feature into shapefile layer. I got the error:

"Attempt to write non-polygon  geometry to polygon type shapefile". Please see my image.

How can I add linestring into layer (shapefile) ?


texts: Group by index column, compare area and add values to new column in Python?

texts:
I have these columns:

index, area_of_poly, cad_admin

I have to group by index (it is a normal column called index) in order to take the rows that have the same value.

#all the ones, all the twos, etc

Some of them (rows) are unique though.
About the ones that are not unique now:

What I have done so far:
I have to check with a group by which of the groups have the largest area and give its respected cad_admin value to the others in its group in a new column called cad_admin1.
The unique values are going to still have the same value they had in cad_admin in the now cad_admin1 column.

more info:

the dataset: http://www.mediafire.com/file/x4q5k7xuztq6o3w/p.zip

import geopandas as gpd
inte=gpd.read_file('in.shp')

inte['index'].value_counts()[inte['index'].value_counts()>1]

359    9
391    8
376    7
374    6
354    5
446    4
403    4
348    4
422    4
424    4
451    4
364    3
315    3
100    3
245    3

inte["rank_gr"] = inte.groupby("index")["area_of_poly"].rank(ascending = False, method =
"first")
inte["key1_temp"] = inte.apply(lambda row: str(row[""]) if row["rank_gr"] == 1.0
else "", axis = 1)
inte["CAD_ADMIN_FINAL"] = inte.groupby("index")["key1_temp"].transform("sum")
print (inte[["area_of_poly", "index", "CAD_ADMIN", "CAD_ADMIN_FINAL"]])

This code as you will see produces some errors like:

TypeError: 'str' object cannot be interpreted as an integer

During handling of the above exception, another exception occurred:

KeyError: ('', 'occurred at index 0')


texts: Understanding target layer and source layer in Select By Location of ArcGIS for Desktop?

texts:
I have a lab assignment that i would need some help with. I have a 25km buffer that's in a circle. Within that i am supposed to make a map were the berries that people can pick are. The assignment is to put in some "demand" that needs to be fullfilled to be in the buffer. I have a layer with vegetation that contains many different vegetation codes. I have made a Selection > Select By Attributes and exported the data with only the vegetation codes that i need.

The problem comes when i run Selection > Select by location. What is the target layer and what is the source layer? I don't understand what the difference is? What is the "acive" layer and what is the "search layer"?

EDIT:
Okay i will try and explain. I tried to upload pictures to be even more clear but i am not allowed to post links apparently. This would be SO MUCH easier to explain if i could post pictures.

First i add my layers of Sweden, the city, islands, roads etc.
I create a pointlayer shapefile to buffer from called R25 and edit in a point right at the city core.
Then i go to geoprocessing > Buffer and chose to buffer R25 to 25km. I save the bufferzone as Buffer25
Then i go to select features and click "Luleå kommun" and then the contour of the city area turns blue. I export this as Data > Export data > Lulea
I go to analysis tools > Extract and chooses clip. There i chose Buffer25 as indata and Buffer25_Lulea as out data. Now there is a new shapefile called Buffer25_Lulea
Now there is a layer called Buffer25_Lulea with a radius of 25km from city core with all the areas that doesn't belong to the city cas been cut off.
Now i have to add the layer called "vegetation". The vegetation types im interested in is vegcode 41, 51 and 52. So i go to Selection > Select by attributes and write "VEGCODE" = 41 OR etc.I save this as Data > Export data and saves is as veg_41_51_52.shp
Now i go to Selection - Select by Location and the method intersect. HERE IS THE PROBLEM!! What layer is the target file and what layer is the source file????? I want all the berries that are in the city to be shown as "Berries Lulea"


texts: Do you need to filter a final GPS [NOT PSEUDORANGE] data with Kalman Filter again?

texts:
I have a projected data in northing, easting, elevation that came from a GPS receiver fitted on a Pipeline Inspection Gauge (PIG). (I do not have the raw data at all)
This PIG collected the GPS data by going through an old pipeline which is buried a few meters under the ground at most places. (There are no drawings showing the pipeline original alignment)
The northing, easting, elevation data looks fine, i.e. it plots as a nice line without any oscillations or zigzag pattern.
I have have been told to correct this data, i.e. to make sure that the alignment we see is the actual alignment of the pipeline underground.
I do not have any additional control point or Inertial Measurement Unit (IMU) data, i.e. the bends depicted by the GPS data may not be a true bend at all.
I have been told that you can use the Kalman filter to filter the northing, easting, elevation data.
So, my question WHY? Why do you want to filter the GPS data that has already been filtered inside the receiver (and probably with a kalman filter)?
Moreover, I have been arguing that we need additional information to check the accuracy of the data and the northing, easting, and elevation information ALONE cannot be used as a input for further filtering.


texts: Add column to point data of whether point falls inside polygon

texts:
This seems like it should be fairly simple, but have been having quite a bit of trouble. I am no expert in ArcGIS for Desktop

I have:
one layer of points, and 4 polygons, which partly overlap

I can select points based on a polygon, however what I want to do is:

I would like to add 4 columns to the points attributes table (to export), with 1s and 0s of whether the point falls inside or outside each of the polygons separately. Is that possible?


texts: Using ModelBuilder to iterate through LAS files using custom LASToolbox?

texts:
I'm trying to use ModelBuilder to iterate through a folder containing 25 usgs .las files and run 4 custom geoprocessing tools from the LAStools toolbox available on github.  I am able to successfully run each tool individually outside of ModelBuilder.  However, when I attempt to put the tools into a ModelBuilder workflow, the ModelBuilder says 'completed' without executing any of the tools.

Below is a screenshot for the workflow with one of the tools which is where I'm starting for now.  I included the tools parameter dialog box which I am leaving default (I used the same parameters when running on its own successfully).

Any ideas on what could be causing this?

I've tried to specify an output name for each tool and still did not execute.  I need to use this tools output as input for the next tool and so on.

FYI, first time using iteration tool in ModelBuilder.


texts: QGIS Print Composer - SVG output displays objects outside the frame

texts:
On Windows 7. QGIS 2.6.0 from OSGEO network installer 64bit (and checked with latest dev version)

I am preparing a print composer map displaying only vector data (shp). Data in the map canvas and the print composer are of same epsg 32624 projection. Input shp data is c. 1-2 mb. The map in the print composer is displaying a subset of the map data. In item properties of the map item, I have pressed the "View extent in map canvas" button to assure, that the map canvas and map item is of roughly the same extent. I've given the map item a thick frame to distinguish it clearly.

The SVG output is 76 kb of size, so far from huge. Nevertheless, the output SVG displays vector map data outside of the map item frame. The problem does not occur in image (png etc.) or pdf output, and neither using raster input data, eg. GeoTIFF.

I've searched the qgis tag, but haven't found a completely similar issue reported, as e.g. memory doesn't seem to be a problem as the input or output file sizes a both very small.

Karl


texts:
When I want to load sld file in geoserver appears this error:

javax.net.ssl.SSLHandshakeException: sun.security.validator.ValidatorException: PKIX path building failed: sun.security.provider.certpath.SunCertPathBuilderException: unable to find valid certification path to requested target
sun.security.validator.ValidatorException: PKIX path building failed: sun.security.provider.certpath.SunCertPathBuilderException: unable to find valid certification path to requested target
PKIX path building failed: sun.security.provider.certpath.SunCertPathBuilderException: unable to find valid certification path to requested target
unable to find valid certification path to requested target

I saw that is something about https.

This is the url:

https://geo.gob.bo/geoserver/ows?SERVICE=WMS&TRANSPARENT=TRUE&VERSION=1.3.0&SLD_VERSION=1.0.0&REQUEST=GetLegendGraphic&CRS=EPSG:3857&SLD=https://geo.gob.bo/mapfishapp/ws/sld/geodoc439edd21e7b1c8ea5a24f35116aef572.sld&LAYER=AGETIC:estado_entidades2017&FORMAT=image/png&LEGEND_OPTIONS=forceLabels:on;fontAntiAliasing:true&SCALE=10000


texts: IOError: Could not open report template from ArcPy?

texts:
I've edited the original post to replace a snippet of code that was producing the same error I was getting on my original code. Below is the original code that was working and then stopped after someone modified the rlf file. The script could no longer open up the template even though it was closed. The next day it ran with no issue. Does anyone know if a memory lock causes this error even if the rlf file is closed?

# Import modules

import arcpy
import os

# Create program variables

# workspace
mxd = arcpy.mapping.MapDocument(r"C:\RS_Data\workspace\BATCH_Carts_SF\BATCH_Carts_SF.mxd")

try:
    # dataframe
    df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]

    # field for SearchCursor
    numField = 'Route'

    # feature layer for cursor
    lyr = arcpy.mapping.Layer(r"C:\RS_Data\Workspace\BATCH_Carts_SF\gisdb\layers.gdb\Service_Location")

    # create empty set to hold unique route values
    rtVals = set()

    # Create searchCursor
    cursor = arcpy.SearchCursor(lyr)
    for row in cursor:
        rtVals.add(row.getValue(numField))

    # create variable to get max value for output name

    scenario = int(max(rtVals))

    path = "S:\\Shared\\RouteSmart\\Automated BIC and Carts\\Test\\Carts\\"

    outFolder = os.path.join(path,("{0}{1}{2}".format(path, scenario, "_Truck\\")))

    outName = os.path.join(outFolder, ("{0}{1}{2}".format(outFolder, scenario,"_Truck_Carts_Report.pdf")))
    print outName

    rlf =r"N:\Common\Report Templates\CO_SF_Carts.rlf"
    arcpy.mapping.ExportReport(lyr,rlf,outName, "USE_RLF")

except arcpy.ExecuteError:
    msgs = arcpy.GetMessages(3)
    arcpy.AddError(msgs)
    print msgs

print "done"

And here's the error from a snippet of the export report only part of the code (sorry I can't get the original error now that the script is working):

Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 323, in RunScript
    debugger.run(codeObject, __main__.__dict__, start_stepping=0)
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\__init__.py", line 60, in run
    _GetCurrentDebugger().run(cmd, globals,locals, start_stepping)
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\debugger.py", line 654, in run
    exec cmd in globals, locals
  File "N:\Common\Script\Drivers\export_report_test.py", line 1, in <module>
    import arcpy
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\utils.py", line 181, in fn_
    return fn(*args, **kw)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\mapping.py", line 515, in ExportReport
    return report_source._arc_object.ExportReport(*gp_fixargs((report_layout_file, output_file, dataset_option, report_title, starting_page_number, page_range, report_definition_query, extent, field_map), True))
IOError: Could not open report template

This is the additional error message I get when using try and except:

PYTHON ERRORS:
Traceback info:
  File "N:\Common\Script\Drivers\CO_051_Export_Carts_Report.py", line 5, in <module>
    import arcpy

Error Info:
<type 'str'>

ArcPy ERRORS:

ArcPy ERRORS:


texts: Can OpenLayers layers in QGIS have transparency?

texts:
I often would like to set the basemap back by lowering the transparency. I would love to have this option in QGIS for OpenLayers layers, such as the stamen layers, google, etc.

Is this possible?

I don't see the ability in the layer properties or within the raster transparency plugin.


texts: What is difference between FID and OID in ArcGIS platform?

texts: What is the main difference between OID and FID in Esri software?


texts: Using Google Fusion Tables to read Wordpress database... which in turn creates map. Can it be done?

texts:
Is it possible for Google Fusion Tables to read my Wordpress database? If so, how?

I have custom wordpress fields which contain addresses for several hundred companies. Ideally I want Google Fusion Tables to read the database.

Hopefully this will overcome two major problems 1. Avoid the geocode limit (currently done using a time delay, which isn't ideal). 2. Obfuscate my data so it's protected.

Really need someone's help with this.


texts: Batch Convert GRID to TIFF with colormap and RGB

texts:
There are lots of questions on here regarding raster conversion, but I can't find anything specific to my needs (involving output as RGB and using a Colormap).

I have a folder full of ESRI GRID format files that I need to convert to TIFF. The part that I'm finding difficult is having the output come out as RGB (I hope to eventually convert these to tiles). When I export one of my grid files as in this first image, it works great. If I just try to straight convert the GRID to tif using something like

arcpy.RasterToOtherFormat_conversion

it creates an image that looks all black to irfanview or any other image viewing program (possibly because the GRID files are binary 1 or 0 ?). It needs to show values of 1 as a color (best if I can make it always the same color of my choosing) and 0 as no color (if possible).

This works great when I do it this way ...

But I want to do this programmatically because there are quite a few of them. Does anyone know of a way to achieve the results I get from Data Export using arcpy and Python (maybe with Environment variables). What am I missing ?

I've noticed that if I have a clr color map file named after the input grid file, and check the "Use Colormap" box, I can get consistent output. This is the color map I was testing with.

0 0 0 0
1 56 168 0

UPDATE: I have tried some of very helpful the suggestions provided in comments below by Michael Miles Stimson and can get the desired result using CopyRaster_management with a colormap file (the input can be 1 or 3 bands, it is the output that needs this file to inform the output). I had some trouble with NoData values (in a separate question) but Chris W helped me with that. I tried using gdal_translate -Expand RGB, but I didn't see any way to specify the exact colors I wanted in the output.


texts: Clearing Calculate Value before Model Run in ArcGIS ModelBuilder?

texts:
I am using ArcGIS 10.2 for Desktop.

I have a model which contains several branches, which branch should run is based on the results of an earlier stage in the model. If the first stage of the model returns a Null result Branch A should run otherwise Branch B should run.
I check for a Null result by using the Get Raster Properties tool ALLNODATA type.

At the start of each branch I created a Calculate Value tool which checks the Get Raster Properties and returns a 1 or 0 based on if that branch should run, and then uses the output as a precondition for the rest of the branch.

This sort of works but once the tool has been run once the various outputs are pre-set and instead of running the "Get Raster Properties" and "Calculate Value" tools again ModelBuilder jumps over to run the branches which meet their preconditions.

Is there a way to set all the results of Calculate Value to 0 before each run of the model so that branches are not run until the Calculate Value tool has run?


texts:
I have the MBtiles but I want to turn the contours within that file into a PDF image.

is this possible or am I wasting my time?


texts: Creating polyline which is average of multiple other polylines?

texts: I am trying to use ArcGIS for Desktop to create an "average" polyline from several lines. For example if you have multiple hurricane paths and would like to plot the average position of all of the paths. I imagine you could do this if you got the y coordinates for specific x coordinates and worked out the mean of the y coordinates. I am not good at scripting and was wondering of there was a tool that could do this.


texts: Update QGIS didn't solve error "NUMPY driver was compiled against GDAL 2.1"

texts:
When I tried "zonal statistics" with QGIS, I encountered the following error.
ERROR 1: NUMPY driver was compiled against GDAL 2.1 but current library version is 1.11

I could find people who had same problem,

QGIS 2.18: Error in NUMPY and GDAL 2.1 [duplicate]
QGIS 2.18.5 update breaks on Ubuntu 16.04 [closed]

so, I updated QGIS 2.18.5 to the latest version.
But I still have the problem.

Here's what I did.

Uninstalled QGIS 2.18.5 by "sudo apt-get remove qgis"
re-installed the latest version, that is, 2.18.7 as the follows.

sudo apt-get update
sudo apt-get -f install
sudo apt-get install qgis python-qgis qgis-provider-grass

Re-installed QGIS says its version is "2.18.7".
I think the above problem should have been solved already,
but I got exactly the same error.

Any ideas ?

--- Update ---

It seems that "sudo apt-get qgis" loads gdal-1.11 even though gdal-2.1 is installed. How can I make it load gdal-2.1 ?


texts: Error importing XML workspace document to file geodatabase

texts:
Ok so I have a large SDE geodatabase that I exported to an XML Workspace Document. I tried to importing it again into a file geodatabase to test if the export worked correctly and got the following error:

Import xml data failed
FDO error: -2147186428
[XML parse error: Illegal xml character
[Line: 23612914, column: 1043].]

Could this be from an invalid character somewhere in my data? Is there a list of invalid characters for XML workspace documents? And is there an easy way to get to the line and column in the error message to see what it is?

Also are there any file size limitations, my XML is nearly 12gb? Or could it be possible there was just an error in the exporting process?


texts: GetParameterAsText() error in python script

texts:
i have a huge problem with an error . when i use getparameterastext and use it for searchcursor , i encounter to an error .for example the below code

import arcpy

featureClass = arcpy.GetParameterAsText(0)

populationField = arcpy.GetParameterAsText(1)

rows = arcpy.SearchCursor(featureClass)

the error is :

  Runtime error : ERROR 999999: Error
  executing function.

the error is different in pythonwin . the error is :

  Traceback (most recent call last):   File "", line
  1, in    File "C:\Program
  Files\ArcGIS\Desktop10.0\arcpy\arcpy__init__.py", line 804, in
  SearchCursor
      return gp.searchCursor(*args)   File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing_base.py", line
  357, in searchCursor
      self._gp.SearchCursor(*gp_fixargs(args))) RuntimeError: ERROR 999999: Error executing function.


texts:
I have two polygons where one should cover the other.  They were created in qgis using snapping, however when I checked their topology with postgis it turned out one point on the inner polygon registered as outside the covering polygon.  This seems to regularly happen when I create an inner polygon with an edge that lies on an outer polygon.

When I use ST_Equals on the point supposedly outside the polygon and the closest point on the outer polygon it returns true, but when I check whether the polygon and the point are touching it returns false.

Here's the query I'm running:

--CTE with points that are outside the outer polygon (only one in this specific case)
WITH outsidepoints AS (
SELECT dump.geom
FROM (
    SELECT geom,
           ST_Covers((SELECT geom from zone WHERE gid = 403), geom)
    FROM ST_DumpPoints((SELECT geom from zone WHERE gid = 416))
    ) AS dump
WHERE st_covers = false
)

--Check if closest point on outer polygon and "outside" point on inner polygon are equal
--Check if "outside" point touches outer polygon
SELECT ST_Equals(geom, ST_ClosestPoint((SELECT geom from zone WHERE gid = 403), geom)),
       ST_Touches(geom, (SELECT geom from zone WHERE gid = 403))
FROM outsidepoints

The result is T, F

Is this some kind of floating point error?  Is there a tolerance parameter I need to set for these calculations?


texts: How to generate jVectorMap for United Kingdom (UK)?

texts:
I have shapefile for UK , i need to convert this to the jVectorMap
Here is my script

python ^
   converter.py ^
   ./shapefile/map.shp ^
   ./output/jquery-uk.js ^
   --width 400 ^
   --where "ISO = 'GBR'" ^
   --language en ^
   --country_name_index ? ^
   --country_code_index ? ^
   --minimal_area 40000 ^
   --buffer_distance -0.5 ^
   --simplify_tolerance 10000 ^
   --longitude0 ? ^
   --name UK

how to find the values(?) country_name_index, country_code_index and longitude0?

I have referred the following pages,

https://stackoverflow.com/questions/11068645/how-to-generate-a-new-map-for-jvectormap-jquery-plugin

http://jvectormap.com/documentation/gis-converter/


texts: Import list of addresses into MapPoint for route planning

texts:
I have a list of 40+ points for each work day in a month.
I need to somehow calculate the distance between each of those points for each day.

I have a list of addresses in a database, house number, street, city, zip and state are all in separate columns so I can generate any output using some simple programming.

But I am stuck at importing addresses to MapPoint. I have tried import wizard but all it does is plots points on map and does not add to route planning.

Maybe there is an alternative to this? Or other way around it?


texts: Restricting access to GeoServer webapp to Localhost in Tomcat?

texts:
I am currently running Tomcat 7 on a Windows Server 2008 machine.  I have GeoServer 4.2.1 installed as a war.  What I need to do is restrict access to GeoServer Manager only to Localhost.  As it currently stands you can access the Manager anywhere by going to

https://www.mydomain.com/geoserver

I would like to only give access to Geoserver Manager by means of

https://127.0.0.1/geoserver

I have updated by context.xml file to limit access to Tomcat Manager via localhost with this:

<Valve className="org.apache.catalina.valves.RemoteAddrValve"allow="127.0.0.1" deny=""/>

Would I do something similar?


texts: Converting 2D shapefile to Equal Area Projection?

texts:
I have thousands of shapefiles with coordinates in lat/long (-180 to +180 x, -90 to +90 y), that contain polygons of which I wish to calculate the area.

I'm using the Python GEOS/OGR/GDAL API.  If I call the Area() function for each geometry, I'll get the 2D area, or the area on a Mercator projection (basically).  Since the polygons are really on a sphere, I want to convert them to an equal area projection before I calculate the area.  I was wondering if anyone had any experience with this on here.

Some reading of the documentation suggested there were functions somewhere in OSRSpatialReference to transform a point or geometry from one projection to another, but it wasn't quite clear to me exactly how to use this or what projection to use.


texts: managing postgresql database table using the pgadminIII sql tool

texts:
i have created a database using postgresql 9.1 and a table containing parcel information for a town council  database on a windows 7 OS. i created the geometry column using the SELECT AddGeomatryColumn function and left the geometry column empty at first. now i want to insert the coordinates of each parcel using the sql tool from the general plan which is in a hardcopy. which query must i use to insert the coordinates one by one? i think i can insert coordinate for each parcel by specifying the erf number of the parcel.

any help please


texts: Identifying reference systems from coordinates in WKT string?


texts: Changing symbols not found for selection in ArcGIS web map?

texts:
Here are two screenshots which will show what I am expecting to get. I am looking for change symbols in the more options tile enlisted in the content tab.

I am expecting to be led to another window which may look like this

Can you advise me in this situation?


texts: Choropleth Mapping Using Quantile Method in CartoDB

texts:
My understanding of choropleth maps that use the quantile method is that they will assign an equal number of items to each bucket (see http://www.ncgia.ucsb.edu/cctp/units/unit47/html/comp_class.html).

A recent map I made with CartoDB generated a map that did not do this. Instead, most items were assigned to one bucket.

Is CartoDB using a different type of quantile method than what I described?


texts: Calculating interior angles of polygons or lines in QGIS

texts:
The problem
In QGIS, it's easy to calculate areas of polygons or lengths of lines using QGIS expressions with area ($geometry) / $area or length($geoemtry) / $length respectively. However, there is no expression to measure angles based on the three points that define it.
The question
Given a polygon of any shape in QGIS, how can interior angle be measured apart from the manual measuring tool (toolbar icon)?
What I tried
What sounds quite easy at the first glance is more complicated. QGIS expressions offer a few angle-functions as:

line_interpolate_angle
angle_at_vertex
main_angle
azimuth

However, there is no get_angle_at_vertex(geometry,index) function so that you could calculate the angle of any polygon- or line-geometry at a certain vertex and neither an angle (point1,point2, point3) function.
Given the points 1, 2 and 3: how to get the angle at point 2, formed by the line connecting these three points on the left side in drawing-order?


texts:
I ran into a file today that is in a format I'm unfamiliar with, and I'm hoping that someone can identify it.

The file extension is CTZ, but the only file type I found on Google with that extension was the "Custody Toolbox Zip" file extension (which is software unrelated to GIS).

Here are a couple of lines from the file (all of the lines are in this format). The whitespace is all space characters, not tabs.

00000 000       0000CANADA                          +1285448500+4308048100+0001891000+0728293500+7050177200+0002341100+0557155000+7257869800+0001302700
10000 000       0000NEWFOUNDLAND - TERRE-NEUVE      +0024419500+0340743500+0001395400+0014446000+0243530500+0001685800+0009974000+0097218800+0000974700
10000 000       6999NFLD./T.-N. CT - SR S/TOTAL     +0007789500+0138134200+0001773300+0004317000+0093599600+0002168200+0003472500+0044536500+0001282500

...

59935A9350005.002412VICTORIA                        +0000188000+0003855300+0002050700+0000096000+0002484500+0002588000+0000092000+0001372500+0001491900
59935A9350006.002413VICTORIA                        +0000109500+0002311900+0002111300+0000056000+0001599600+0002856400+0000053500+0000709900+0001326900

Thanks!

Update: After some more Googling, I found this mailing list post of a person with very similar data, but I couldn't garner the format from that conversation.

Update 2: The data was acquired via the Statistics Canada Data Liberation Initiative. It sounds like it came from an FTP site that requires authentication, so I can't provide access to the files.


texts: Population data at OpenRouteService QGIS Plugin

texts: Does anyone have any experience of how reliable the accessibility analysis and the number of inhabitants in the catchment areas are for German cities?


texts: ArcPy / ArcMap 10.2: How to report data to user?

texts:
I'm developing an arcpy addin with ArcMap 10.2 that outputs a feature layer based on user input.

How can I display attributes from this layer in a MessageBox or other reporting mechanism?

I don't want to enable data driven pages, so dynamic tables don't seem to be an option.


texts: I need to configure geoserver proxyHost even standing on my local machine?

texts:
I'm with geoserver on my local machine running on port 8080 on jetty which is standard in windows installer.
I want to use my GetFeatureInfo WMS layers, but I have the following questions:

I need to configure proxyHost even being on my local machine?

Can I use ProxyHosto GeoServer running with Jetty as the default installation for WindosEstou with geoserver on my local machine running on port 8080 on jetty which is standard in windows installer.
I want to use my GetFeatureInfo WMS layers, but I have the following questions:

I need to configure proxyHost even being on my local machine?

Can I use ProxyHosto with GeoServer running on Jetty as is standard for Windows installation or do I have to install it in Tomcat?

This is my code, when I click is carrying too much time and nothing appears:

** EDITED CODE **

<!DOCTYPE html>
<html lang='en'>
<head>
 <meta charset='utf-8' />
 <title>My OpenLayers Map</title>
 <script type='text/javascript' src='OpenLayers.js'></script>
 <script type='text/javascript'>

 OpenLayers.ProxyHost = "proxy.cgi?url=";

 var map;

 function init() {

    var bounds = new OpenLayers.Bounds(
                -51.280762, -30.241074,
                -51.020828, -29.969619
            );
            var options = {
    controls :[new OpenLayers.Control.Navigation(),
            new OpenLayers.Control.PanZoom()],
                maxExtent: bounds,
                maxResolution: 0.00106037109375,
                projection: "EPSG:4326",
                units: 'degrees'
            };
    map = new OpenLayers.Map('map_element', options);

var wms = new OpenLayers.Layer.WMS(
            "Paradas de ônibus",
            "http://localhost:8080/geoserver/wms?",
            {
               layers: 'Teste:Stops',
               isBaseLayer: false,
               styles : '',
               format :"image/jpeg"
            }
        );

    map.addLayers([wms]);

    //GetFeatureInfo Simples
    info = new OpenLayers.Control.WMSGetFeatureInfo({
        url: 'http://localhost:8080/geoserver/wms?',
        title: 'Identify features by clicking',
        queryVisible: true,
        eventListeners: {
            getfeatureinfo: function(event) {
                map.addPopup(new OpenLayers.Popup.FramedCloud(
                    "chicken",
                    map.getLonLatFromPixel(event.xy),
                    null,
                    event.text,
                    null,
                    true
                ));
            }
        }
    });
    map.addControl(info);
    info.activate();

    var point = new OpenLayers.LonLat(-51.22,-30.08);
    point.transform(new OpenLayers.Projection("EPSG:4326"),
    map.getProjectionObject());
    map.setCenter(point, 3);

    map.addControl(new OpenLayers.Control.LayerSwitcher({}));

 if(!map.getCenter()){
    map.zoomToMaxExtent();
 }
}

 </script>

</script>
</head>

<body onload='init();'> <!-- Chama a função js init() -->

<!--Elemento HTML onde o mapa é exibido -->
 <div id='map_element' style='width: 800px; height: 800px;'>
 </div>

</body>
</html>


texts: Snapping in Esri javascript api does not work on mobile device

texts:
I have an Esri map loaded using the javascript api. I try to enable snapping, where a point should always snap to a given line.

I have a snap point symbol that will trace the line to snap to when I drag the edited point around. When I release the mouse, the edited point should snap to the line.

This works perfect on desktop, but when I try on iPhone 5 and HTC One, the snapping is not activated. The snap symbol is never shown when I drag the point, and it is not snapped when I release.

This is the code for how I enable the snapping:

map.enableSnapping({
    snapPointSymbol: snapSymbol,
    tolerance: 20000,
    alwaysSnap: true,
    layerInfos: [{ layer: map.graphics, snapToPoint: true, snapToVertex: true, snapToEdge: false}]
  });

Emulation testing:

When I try to emulate this in Chrome developer tools, I see the same problem for iOS devices, but not for Android or any other types of device. It seems like the ESRI code is sniffing the user agent string, and if it finds an iOS specific value in the user agent, the snapping is disabled.

Real device testing:

Reproduced on iPhone 5 (iOS 7) and HTC One.

How can I get snapping to work on mobile device?


texts: Using Collector with ArcGIS for Personal Use

texts:
If I buy a personal use account will I be able to use collector to collect data in the field?

For context, I would like to be able to collect data points in the field to make my own hiking and trail maps.

If I will not be able to use collector, what would the community suggest?


texts: NDVI chart using sentinel 2 image in Google Earth Engine

texts:
I want to make a chart in google earth engine based on monthly NDVI value of one plot. For that i am using the following codes

  function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = ee.Number(2).pow(10).int();
  var cirrusBitMask = ee.Number(2).pow(11).int();

  // Both flags should be set to zero, indicating clear conditions
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
             qa.bitwiseAnd(cirrusBitMask).eq(0));

  // Return the masked and scaled data.
  return image.updateMask(mask).divide(10000);
}

// Plot Landsat 8 band value means in a section of San Francisco and
// demonstrate interactive charts.

//var sanFrancisco = ee.Geometry.Rectangle(-122.45, 37.74, -122.4, 37.8);

var boundary = ee.FeatureCollection('ft:1T-XrxMki2k49fhQdhMxesWKBoyyUJiKN70pErh9z');// Boundary
var gtpoints = ee.FeatureCollection('ft:1ZvFeFO5jHGo3VExlKNc-kdL-3hmJSeyfLzqLvHhy');//Co-ordinates
var plots = ee.FeatureCollection('ft:1pZJ7v8ek4Z1GP7MxZNA3EhtRVGs5MDgEP6ngwc1v');//Plots
var startdate, enddate;
var NDVI = [];
var ndvi = [];
var month = [];
var mean;
for(var i = 1; i<=12; i++)
{
  startdate = '2017-'+i+'-01'
  enddate = '2017-'+i+'-28'
 // print(startdate, enddate);

var image = ee.ImageCollection(sent2image
.filterDate(startdate,enddate)
.filterBounds(boundary)
//.map(maskS2clouds)
.sort("CLOUD_COVERAGE_ASSESSMENT")
.median());

var mosaic = image.mosaic();
var clip = mosaic.clip(field);

var ND = clip.normalizedDifference(['B8','B4']).rename('nd');

var mean = ND.reduceRegion({
  geometry: field,
  reducer: ee.Reducer.mean(),
  scale: 10
});

mean = mean.get('nd')

month.push(i)
NDVI.push(mean)
}
print(NDVI)
print(month)

var chart = ui.Chart.array.values(ee.Array(NDVI),1)
print(chart)

I am perfectly able to get the mean ndvi values in an array. But i am not able to produce the chart for the same
Sharing the link also for my code: https://code.earthengine.google.co.in/00fd8691e53c8ea0c7968276e5f4f5ac


texts: Geoserver MVT - How to serve 'id' attribute per vector tile spec?

texts:
I am trying to work with MVT tiles served by a geoserver instance.

The vector tile spec contains obligatory fields for type and geometry, and optional fields for id and tags.

The MVT tiles as served by the GeoServer instance leaves the id field as undefined, which isn't ideal when working with hover / selection workflows.

Is there a way to configure the GeoServer instance to automatically generate a unique id for populating the id field?

Example feature data as returned from mapbox queryRenderedFeature:

{
  type: "Feature",
  properties: {
    col_id: 'uniqueId',  <-- the table data contains a unique id
    col_1: 'some data',
    col_2: 'more data'
  },
  id: undefined,  <-- this "outer" MVT id is undefined
  geometry: ...
}


texts:
I have made a raster elevation map of Greenland, and have information of the region, cities and more as vectors. The information is not in the same place as my elevation map.

I have selected EPSG:5938 on the bottom right in QGIS and on my informations and elevation map I have set the CRS to the same like the picture

How do I get it all on the same position?


texts: Geoserver WMS layer becomes pixelated in OpenLayers

texts:
I added a Geoserver WMS to OpenLayers, and if I zoom, the layers become pixelated. I've created the layers and SLDs (for symbology)  with QGIS. Is the problem with the zoomlevels or elsewhere?

var map = new OpenLayers.Map('terkep_helye',
   {
     projection: "EPSG:4612", // EPSG:3857 = Spherical Mercator
     maxExtent: new OpenLayers.Bounds(-179.76129999999998,-60.60800000099994,179.8569000000001,59.38540000100005),
   });
 // hozáadunk egy OSM réteget
 map.addLayer(new OpenLayers.Layer.OSM());
  // Google hybrid
 map.addLayer(new OpenLayers.Layer.Google("Google Hybrid",{type: google.maps.MapTypeId.HYBRID, numZoomLevels: 12}));
 // Google satellite
 map.addLayer(new OpenLayers.Layer.Google("Google Satellite",{type: google.maps.MapTypeId.SATELLITE, numZoomLevels: 12}))
 // létrehozunk egy WMS réteget var wms = new OpenLayers.Layer.WMS( "Földrengések",
    'http://localhost:8080/geoserver/ujProjekt/wms',
    {
      layers: ['rengesekk'],
      isBaseLayer: false,
      transparent: true,
      srs: "EPSG:4612",
      numZoomLevels: 24
    },
    {

      singleTile: true
    });
var wms2 = new OpenLayers.Layer.WMS( "Szubdukciós vonal",
    'http://localhost:8080/geoserver/ujProjekt/wms',
    {
      layers: ['tvonal'],
      isBaseLayer: false,
      transparent: true,
      srs: "EPSG:4612",
      numZoomLevels: 12
    }{
      //opacity: 0.5,
      singleTile: true
    });
    var wms3 = new OpenLayers.Layer.WMS( "Intenzitás",
    'http://localhost:8080/geoserver/ujProjekt/wms',
    {
      layers: ['intenzitas'],
      isBaseLayer: false,
      transparent: true,
      srs: "EPSG:4612",
      numZoomLevels: 12
    },
    {
      //opacity: 0.5,
      singleTile: true
    });
    var wms4 = new OpenLayers.Layer.WMS( "Maximális Talajgyorsulás",
    'http://localhost:8080/geoserver/ujProjekt/wms',
    {
      layers: ['maxtalajgyors'],
      isBaseLayer: false,
      transparent: true,
      srs: "EPSG:4612",
      numZoomLevels: 12
    },
    {
      //opacity: 0.5,
      singleTile: true
    });
    var wms5 = new OpenLayers.Layer.WMS( "Domborzat",
    'http://localhost:8080/geoserver/ujProjekt/wms',
    {
      layers: ['domborzat'],
      isBaseLayer: false,
      transparent: true,
      srs: "EPSG:4612",
      numZoomLevels: 12
    },
    {
      //opacity: 0.5,
      singleTile: true
 });
 // hozzáadjuk a réteget a térképhez
 map.addLayer(wms);
 map.addLayer(wms2);
 map.addLayer(wms3);
 map.addLayer(wms4);
 map.addLayer(wms5);
 // a LayerSwitcher kezelőszerv segítsével kapcsolgathajuk a rétegeket

 // kb. Magyarországra nagyítunk
 map.zoomToMaxExtent();
 map.addControl(new OpenLayers.Control.LayerSwitcher);
 map.addControl(new OpenLayers.Control.OverviewMap());

 map.addControl(new OpenLayers.Control.Scale());
 map.addControl(new OpenLayers.Control.Navigation());
 map.addControl(new OpenLayers.Control.ScaleLine({
    maxWidth: 200,
    bottomOutUnits: '',
    bottomInUnits: ''
  }));

}

OpenLayers - térkép több réteggel

Az alap az OSM (OpenStreetMap), rajta egy Geoserver wms réteg


texts: Why does Fetch Python Plugins responded with only the installed ones?

texts: I use QGIS 1.8. Whenever I demand to install plugins through Fetch Python Plugins it displays only those the already installed. I couldn't find for example; point sampling, group stats and photo2shape plugins. Can one help me?


texts: ArcMap Symbology: Assign colors based on a column value in attribute table

texts:
I am very new to ArcMap and I am still learning the basics and terms. If I am not explaining things enough please ask me.

This is what I am trying to do: I have a soils dataset and the Soil map units column I am interested in has values like this
'HfD2',
'FaE2',
'SaA',
'SyB' and so on.

I would like to know if it is possible to assign colors based on the range A - E. i.e., if the map units with 'A' all fall into one color category, 'B' into another and so on. The range is A through E, so I know there's only going to be 5 categories.

I am thinking this should be done with a SQL expression but can not figure out how to do the same.


texts:
I have a single layer of points I need to split into separate layers based on the unique ID. Below is an example of how the points are labeled. First digit is the feature number and after the decimal is a point along that feature. I'd like to separate them based on the feature number (e.g. all points for Feature 1 (1.1,1.2,1.3) are on a separate layer).

1.1
1.2
1.3
2.1
2.2
2.3
2.4
2.5
3.1
3.2


texts:
I am using openlayers to display 3 layers, a google base layer, and 2 layers that i call from Geoserver. The projection i am using is EPSG:3857.
My problem is ......
My overview map location does not tally with the location on my main map. So for instance if i scroll on my main map to north Africa, my overview map shows Australia. The overview map moves when the main map moves but the location is completely different. I tried all sorts of code but i cannot solve the problem.
Does anyone know how i can solve my problem?
One of my codes is below....

map.addControl(
    new OpenLayers.Control.OverviewMap({
        autoPan:true,
        maxExtent: new OpenLayers.Bounds(10264212.20144042, -1368478.0758595483, 17758023.6939145, 3318388.9003530084),
        restrictExtent: new OpenLayers.Bounds(10264212.20144042, -1368478.0758595483, 17758023.6939145, 3318388.9003530084),
        maxResolution: 29272.70114247687,
        //projection: "EPSG:3857",
        projection: new OpenLayers.Projection("EPSG:3857"),
        units: 'm',
        numZoomLevels: 5
    })
);


texts:
I would like to see the geoJSON placemark label once hover it.

I found a few examples, but they refer mostly to providing a pop-up only:

Showing popup on mouse-over, not on click using Leaflet?

https://stackoverflow.com/questions/37408445/leaflet-mouseout-event-on-marker

Show popup on marker hover / mouseover, hide on mouseout, and keep open on click

My code looks as follows:

  job = L.geoJson(data, {

        pointToLayer: function(feature, latlng) {
        feature.properties.myKey = ''+ feature.properties.Owner + '; ' + feature.properties.Address
            return L.circleMarker(latlng, {
            radius:6,
            opacity: .5,
            //color: "#000",
            color:getColor(feature.properties.Contractor),
            fillColor:  getColor(feature.properties.Contractor),
            fillOpacity: 0.8

            });  //.bindTooltip(feature.properties.Name);
        },

            onEachFeature: function (feature, layer) {
                layer._leaflet_id = feature.properties.Owner;
                label = String(feature.properties.Owner)

            function highlightFeature(e){
                var ly = label
                }

                var popupContent = "

" +
                feature.properties.Owner + "
 Number of units: " +
                feature.properties.Units + "
 Client reference: " +
                feature.properties.Client_Ref + "
" +
                '
Local directory
' +
                'Sharepoint
';

                if (feature.properties && feature.properties.popupContent) {
                    popupContent += feature.properties.popupContent;
                }
                    layer.bindPopup(popupContent);
                    layer.on({
            mouseover: function (e) {
                this.setStyle({
                    'fillColor': '#663300',
                    'weight':11,
                });
                highlightFeature(ly);
            },

            mouseout: function (e) {
                this.setStyle({
                    'fillColor': getColor(feature.properties.Contractor),
                    'weight': 2
                });
            }
            /*
            click: function () {
                // TODO Link to page
                alert('Clicked on ' + feature.properties.Title)
            }
            */
        });

            }

            }).addData(data).addTo(map);

where inside the  mouseover: function (e) { I plotted the highligtfeature function highlightFeature(ly);  as added to the onEachFeature: function above.
The map is fine, but I see no effect.
The console says:

  (index):497 Uncaught ReferenceError: ly is not defined
      at NewClass.mouseover ((index):497)
      at NewClass.fire (Events.js:190)
      at NewClass._fireDOMEvent (Map.js:1386)
      at NewClass._handleDOMEvent (Map.js:1343)
      at HTMLDivElement.handler (DomEvent.js:79)

Shall I input the:

     .bindTooltip(label, {permanent: true, direction: "center", className: "my-labels"}).openTooltip();

somewhere in the code too?

Could you help me with it?

I don't want popop opened, only the label.


texts:
I have PostgreSQL 9.2 with Pgrouting extension, but when I want to call a function(pgrouting)
does not work.

Id like to compute the shortest path between two point.

SELECT gid, AsText(the_geom) AS the_geom
        FROM shootingstar_sp('geom_way', 281346, 280088, 0.1, 'length', true, true);

result:

ERROR: function shootingstar_sp(unknown, integer, integer, numeric, unknown, boolean, boolean) does not exist
SQL state: 42883

I did it same way with astar_sp_delta(); I got same result...

Thanks


texts: QGIS 3.0.2 NNJoin doesn't complete. Crashes when close status bar

texts:
I have run a large spatial point to polygon join of 2 million points to 657 polygons. I used NNJoin so that the closest polygon will be assigned to a point if the point does not lie within one.
The process ran for around 6 hours. The status bar showed 100%, but no 'finished' message or result output file was added to QGIS or can be found in my directory.
If I close the status bar on QGIS, I get a Python error "RuntimeError: wrapped C/C++ object of type QProgressBar has been deleted.
Any ideas what the problem is?


texts:
So I am trying to figure out how to understand the Getis Ord results from pysal.

I have a geographically enabled datatframe, and I want to figure out how to find both hot and cold spots from the results.

import pysal

fc = r"C:\temp\towns.shp"
field = "TOTALPOP10"
sdf_poly = SpatialDataFrame.from_featureclass(fc) # polygon

permutationsValue = 999
w_poly = pysal.weights.Queen.from_dataframe(df=sdf_poly, geom_col='SHAPE')
y = sdf_poly[field]
statistics_g = G_Local(y, w_poly,
                     star=True, transform='R', permutations=permutationsValue)

This all works fine and well, now I need to find significant values:

sig = statistics_g.p_sim<0.05
statistics_g.p_sim[sig]

Are these the hots zones?

If so, how do I find the cold spots as well?


texts:
Are there any options to download a GeoPackage from GeoServer through WFS service?

I understand WFS output formats, but I can't find documentation about providing the GeoPackage format.


texts: ArcGIS 10.0 height statistics for polygon features: zonal statistics vs interpolate shape

texts:
I'm performing a fairly routine task, which is to add height information (from a 50m DEM) to survey polygons to work out mean height, as well as max and min.

I've noticed that zonal statistics (Spatial Analyst) is very slow, and depending on the raster analysis settings, can miss out polygons due to raster sampling size. For example, in the image below, the red polygon is missed out of zonal statistics because it's too thin to adequately intersect enough (or indeed any) raster cells.

Interpolate shape -> Add surface information, by contrast, seems much much faster. But I don't understand what's being interpolated here for a polygon feature. Is it just the vertices, or is the interior of the polygon also being compared with my underlying DEM? At what resolution?

Which method would provide more accurate results? And, are there any other solutions. I often find that terrain operations work equally well with rasters containing non-height information, and vice versa.


texts: Expanding existing point grid to cover larger area in ArcGIS Desktop

texts:
I am looking to expand the already existing point grid system to encompass the entire polygon pictured below without changing the existing point locations or label features.

Is there a tool in ArcGIS Desktop that will allow for the expansion of an already existing grid?

I am running on ArcGIS Desktop 10.5.


texts: How do I add a MapProxy overlay to Google Maps API?

texts: I have Mapserver running, reading from a raster file, which I can overlay over Google Maps using its Javascript API. Now, for caching of tiles, I need to use MapProxy URL for Google Maps overlay. How do I do that?


texts: Trying to run a clip_analysis for ArcGIS

texts:
Trying to run a clip_analysis for a buffer analysis I just wrote but I keep getting an TypeError: unsupported operand type(s) for +: 'int' and 'str'

Not sure what I need to do to fix it and be able to run the buffer_analysis first and then the clip_analysis.

Script I have now is

for i in range(len(cities)):
arcpy.Buffer_analysis(sites[i]+'.shp', sites[i]+'_Buffer3000.shp','3000')
print 'site buffered to',i,'buffer 3000'
arcpy.Clip_analysis(openSpace,i+'.shp', i+'_OpenSpace.shp')
print 'open space clipped to',i,'city limits'

The goal is to buffer the sites around the cities feature class and then clip the openSpace to the newly created buffers.

Also both analysis need to be run in the same for loop


texts: Map chopped off across 180degree line from QGIS to Mapbox

texts:
I am currently exporting a geojson data from qgis to mapbox studio. My area map covers and crosses the 180 degree line across the Pacific ocean. I tried exporting using CRS setting EPSG:4326 and EPSG:3832 (Pacific Disaster Center) to no luck as both settings cuts and only shows the first half of the area.
Is there a way that I can export the geojson without being cut off?

Area should be this:


texts: Why can't I see all records in my personal GDB in QGIS?

texts:
I am running QGIS 1.8.  I have a layer from an ESRI Personal GDB.  There are 431 features in the layer.  When you load the layer into QGIS, and look at the Metadata properties, it shows the number of features as 431.

When I open the attribute table, the table only shows 180 features.

This is a layer that had been loaded in straight from the Personal GDB.  No filters or queries had been added to it, so all features should be showing.

I first came across this issue because another coworker was trying to access the same layer using QGIS 1.7.  He was having the same issue.

A couple items of note:

This layer was existing in the PGDB, and then using ArcGIS, I copied and replaced a new version from a File GDB into the PGDB.
I am able to load other layers from the same PGDB into QGIS where all the features are displayed.
The PGDB is a modified MS Access database.  If you open the database with MS Access, then open the table for this layer, it contains all 431 records with attributes.

Here is a screenshot of the layer in QGIS:

Here is a screenshot of the same layer in ArcGIS:

If anyone has some insight into why QGIS is able to see the actual number of features in the layer, but is only displaying a certain number, I would appreciate it.


texts: How many themes in layers box can I create in QGIS 3.0?


texts:
I have created a 1km * 1km grid in QGIS using the MMQGIS plugin and below is a sample of shp@data and the layer

Is there a way to do a spatial join on the grid created?

I assume some grids should have more individual points (lat/long) and some should have no point data within them. I want to be able to identify the "grid_id" and the points present in the polygon such that I can then join the individual point data like below for e.g. in the below-aggregated data, which grid number does the lat and long belong to head(b)

Lat                Long     Num.CT   Num.SuccessCT Num.FailedCT Failure.rate
49.91222         -6.302522      1             1            0       0
49.91251         -6.299382      1             1            0       0
49.91440         -6.315690      1             1            0       0
49.91443         -6.309872      1             1            0       0
49.91501         -6.311369      1             0            1       1
49.91502         -6.310905      1             1            0       0

I tried joining the shapefile with the aggregated data. However, I can't get my head around how to plot and how would leaflet package in R identify the square

Spatial_join_by_location
shapefile join aggregated data

Output post joining the point layer and polygon layer


texts: Google Earth Engine -Iterated image of NDWI

texts:
I computed the NDWI of the whole image collection and I expect to output an image that shows water body whenever it has been flooded , even once , throughout the filtered date.

Somehow the output image only show the result of the last image in the image collection.

This is what I currently have :

// Function to clip image collection
var clipimgcol = function(image){ return image.clip(Border)} ;
// ImageCollection
var l8collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
.filterBounds(Border)
.filterDate('2014-02-05', '2018-02-05')
.sort('system:time_start')
.map(clipimgcol);

// Mask pixels with clouds and cloud shadows
// Surface reflectance products come with a 'pixel_qa' band based on CFMask
// Create function leaving Water and Clear
var maskCloudWater = function(image){
var pixel_qa = image.select('pixel_qa');
return image.updateMask(pixel_qa.eq(322).or(pixel_qa.eq(386)).or(pixel_qa.eq(834))
.or(pixel_qa.eq(898)).or(pixel_qa.eq(1346)) //pixel_qa = Clear
.or(pixel_qa.eq(324)).or(pixel_qa.eq(388)).or(pixel_qa.eq(836))
.or(pixel_qa.eq(900)).or(pixel_qa.eq(1348))); //pixel_qa = Water
};

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

// Compute the Normalized Difference Water Index (NDWI).
var getNDWI = function(img){
  return img.addBands(img.normalizedDifference(['B3','B5']).rename('NDWI'));
};

// map over image collection
var l8ndwi = l8masked.map(getNDWI);
var ndwi = l8ndwi.select('NDWI');
print(ndwi);

I use a threshold of 0.3 for NDWI and tried to iterate the image collection so once the pixel's NDWI is gte 0.3 it returns 1 and the output will be the way I expected, but resulted differently :

var ini = ee.Image(ndwi.first());
function it (current , previousresult){
  var step1 = current.gte(0.3);
  var step2 = ee.Image(previousresult).select('NDWI');
  var laststep = ee.Image.cat([step2 , step1]);

  return laststep;
}
var ndwiit = ee.Image(ndwi.iterate(it , ini)).select('NDWI');
var ndwiViz = {min: -1, max: 1, palette: ['green', 'white','blue']};
Map.addLayer(ndwiit , ndwiViz, 'NDWIIT', false);

I also came up with the idea of remap function, to update the pixel's NDWI value once it is gte 0.3 but i didn't know how to extract band value of every pixel.


texts: gdal_merge.py is not working with NODATA values

texts:
I am trying to merge a couple of GeoTiff files with gdal_merge.py. Example input:

In the input data the nodata value is -32678. Hence I tried:

gdal_merge.py -n -32678 -o DEM.tif -of GTiff na_DEM.tif ca_DEM.tif

and

gdal_merge.py -n -32678 -a_nodata -32678 -o DEM.tif -of GTiff na_DEM.tif ca_DEM.tif

However, for some reason the nodata values are not recognized as the output file  giving me erroneous values where the two files overlap. Below I provided the output from gdalinfo on one of the input files.

Driver: GTiff/GeoTIFF
Files: ca_DEM.tif
       ca_DEM.tif.aux.xml
Size is 14160, 6000
Coordinate System is:
GEOGCS["WGS 84",
  DATUM["WGS_1984",
      SPHEROID["WGS 84",6378137,298.257223563,
          AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
  PRIMEM["Greenwich",0],
  UNIT["degree",0.0174532925199433],
  AUTHORITY["EPSG","4326"]]
Origin = (-119.000000000000000,29.999999999999002)
Pixel Size = (0.004166666666667,-0.004166666666667)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
 Upper Left  (-119.0000000,  30.0000000) (119d 0' 0.00"W, 30d 0' 0.00"N)
Lower Left  (-119.0000000,   5.0000000) (119d 0' 0.00"W,  5d 0' 0.00"N)
Upper Right ( -60.0000000,  30.0000000) ( 60d 0' 0.00"W, 30d 0' 0.00"N)
Lower Right ( -60.0000000,   5.0000000) ( 60d 0' 0.00"W,  5d 0' 0.00"N)
Center      ( -89.5000000,  17.5000000) ( 89d30' 0.00"W, 17d30' 0.00"N)
Band 1 Block=14160x1 Type=Int16, ColorInterp=Gray
  Min=-42.000 Max=4065.000
  Minimum=-42.000, Maximum=4065.000, Mean=836.674, StdDev=815.941
  NoData Value=-32768
  Metadata:
    STATISTICS_MAXIMUM=4065
    STATISTICS_MEAN=836,67402450978
    STATISTICS_MINIMUM=-42
    STATISTICS_STDDEV=815,94103514887


texts: Counts points in polygon for large number of points QGIS

texts:
have a relatively large number of points ~3.5m and want to count the number of these within administrative boundaries (182 polygons). Point data is in a PostGIS db (through a materialized view) polygons are in simple shapefile (not in db)

Initially running the points in polygon just stalled my machine so I've tried breaking down the query using a filter on the points (which gives 150,000 points to count). This still doesn't run. I've tested the shapefile with very small (30) data set and no problems there.

My workround is to save out a shapefile, but i'd like to know why it's not working and if there is an alternate tool in QGIS i can try?  Ultimately I want a basic choropleth map to show how many points in each admin area (polygon)

I


texts: Which version of pgRouting should I use with PostgreSQL 9.1 and PostGIS 2.0 (on win xp)?

texts:
I am new to GIS.
I want to install pgRouting on my system (win-xp 32-bit). I am trying to run routing_core.sql file and got the following error:
ERROR:  incompatible library "C:/Program Files/PostgreSQL/9.1/lib/librouting.dll": version mismatch
DETAIL:  Server is version 9.1, library is version 8.4.
What should I do for proper installing of pgRouting ?
Thanks for your suggestion/comments.


texts: Clipping VRT with mask layer results in too large VRT

texts: i have created with QGIS a VRT holding ~400 raster tiles,  the totaling size of all rasters held in this VRT is around 90GB, the file size of the VRT is 1.9MB. I now want to clip this VRT with a mask layer using the "Clip raster by mask layer" native plugin. As a result i want an VRT again. The operation succeeds with no error message - but the resulting VRT has a file size of 34MB now and it does not load (for ~14h) in QGIS when want to add it to the map canvas. How can i create a VRT with the clipped extent? Is there another, more performant way to achieve this?


texts:
I've got some code which adds menu items with some inactive by default. They only become active when the reverse item (to open the data) above it are triggered.

I am now looking at using a loop to do this rather than having to repeat the code for each change in value.

It's mostly working fine, except for enabling/disabling the menu items.

CODE that is in the loop

                        elif i == 2:
                            a = QAction(row[i], self.iface.mainWindow())# add the last level of menu item
                            a.setObjectName(row[i]) # Assign the name
                            a.triggered.connect(self.loadLayer) # load the layer(s)
                            actions.append(a)
                            submenu.addAction(a)
    #add close option and seperator
                            msg = 'Close '+row[i]
                            cls = QAction(msg, self.iface.mainWindow())
                            cls.setObjectName(msg)
                            cls.triggered.connect(self.unloadLayer) # unload the layer(s`)

self.cls.setEnabled(False) # I tried this but get an error

                            actions.append(cls) #?
                            submenu.addAction(cls) #?
                            submenu.addSeparator()#?

What execatly do these do?

    #end
[![enter image description here][1]][1]

CODE that is not in the loop but works as expected

    menuSewerageNetwork = QMenu('Sewerage Network', self.menuInfrastructure)
    self.actionOpenSewerageActive = QAction('Open Sewerage Active', self.iface.mainWindow())
    self.actionOpenSewerageAll = QAction('Open Sewerage All', self.iface.mainWindow())
    self.actionCloseSewerageActive = QAction('Close Sewerage Active', self.iface.mainWindow())
    self.actionCloseSewerageAll = QAction('Close Sewerage All', self.iface.mainWindow())
    self.actionCloseSewerageActive.setEnabled(False)
    self.actionCloseSewerageAll.setEnabled(False)

    menuSewerageNetwork.addAction(self.actionOpenSewerageActive)
    menuSewerageNetwork.addAction(self.actionOpenSewerageAll)
    menuSewerageNetwork.addSeparator()
    menuSewerageNetwork.addAction(self.actionCloseSewerageActive)
    menuSewerageNetwork.addAction(self.actionCloseSewerageAll)

    self.menuInfrastructure.addMenu(menuSewerageNetwork)

    self.actionOpenSewerageActive.triggered.connect(self.openSewerageAll)
    self.actionCloseSewerageActive.triggered.connect(self.closeSewerageAll)


texts:
I am using QGIS Essen version (2.14.1). I am trying to filter my data based on dates. I am wondering if the "between" operator exists like in SQL?
I'm trying the following query but the results are wrong. Or is it the date format which is a problem?

(
(("all_ws_d_2"  <= '01-09-16') and  ("all_ws_d_3"   >=   '30-09-16'))
or (("all_ws_d_2"  >=  '01-09-16') and ("all_ws_d_3"  <= '30-09-16'))
or (("all_ws_d_2" <=  '01-09-16') and (("all_ws_d_3"   => '01-09-16') and ("all_ws_d_3"   <=  '30-09-16' )))
or (("all_ws_d_2" >= '01-09-16') and  (("all_ws_d_2"  <=   '30-09-16') and ("all_ws_d_3"  >=  '30-09-16')))
)


texts: Align arrow (point layer) to water direction (line feature)

texts:
Assuming that pipes are digitizet correctly (first node as upstream, last as downstream), I need to align the arrows of a point layer to display flowmeter's direction. So far I need to use v.split (GRASS), then populate a new field with this expression (see How to horizontally align irregularly shaped point symbol to vector layer in QGIS?):
(atan((xat(-1)-xat(0))/(yat(-1)-yat(0)))) * 180/3.14159 + (180 *(((yat(-1)- yat(0)) < 0) + (((xat(-1)-xat(0)) < 0 AND (yat(-1) - yat(0)) >0)*2)))

and then this angle can be assigned to nearest point with refFunctions plugin:
geomdistance('pipes', 'Angle', 0.001)

Since the flowmeters will be moved to different pipes, I'm looking for a different procedure that will change automatically the angle, I'm open to advices.
What comes to my mind is:

buffer a circumference (radious like 5 meters) around each point,
select the segment inside buffer
apply the angle formula to this segment and then write the value as point attribute

Questions are: can it be done? Suggestions in how to do it (preferably with field calculator)?


texts:
I am trying to add a couple of HDF rasters into a mosaic dataset, however no matter what I do I keep getting this error:

This is a screenshot of my inputs:

The raster processing and mosaic post-processing default options were left.

Any idea what the problems is?  Is there something that needs to be done with the .H5 files before adding them into the mosaic dataset??


texts: Finding sum of points values within distance from paths

texts:
I'm willing to plot public transport data for a university project, and I'm currently trying to find out roughly how many people live with a reasonable distance from a set of bus routes.

I have my bus routes set up as paths, like the green path below, and population data per small area as an attribute for each pink point.

What I'm trying to do is create a buffer around each bus route of roughly 1 km, select all the pink points within the buffer area, find the sum of the "Population" value of each, and save that in a file. However, I have around 600 routes I'm trying to find this info for, and I haven't been able to find any way to batch process all of my routes.


texts: Problems with import of QGIS generated PDF or SVG files in Adobe Illustrator

texts:
I have some issue with importing PDF or SVG files generated in QGIS 3.16 Hannover (with the option to export everything as vectors, no compression, no optimization for smaller file size, incl. all layers, attributes etc.)
In this particular case, I used a digital terrain model to generate a layer of 19 different contour polygons from 0 to 4500 meters in 250 m steps with different shades of grey. In QGIS 3.16 everything appears to be just as I want it. I use the "new print layout" option to generate a map and exporte to SVG or PDF. If I just open the GeoPDF file in a standard PDF reader it looks like this, which is what I want (it's a random section of the original PDF).

However, importing and opening the same file in Adobe Illustrator gives me something completely unusable:

The same goes for the SVG files, just some unrecognizable garbage.
Any ideas what is going on and how to import a svg/pdf file into Adobe Illustrator with the layer and attributes correctly incorporated? I.e. in this case the 19 different contour polygons with altitude and color info, so that I can continue processing it in AI.


texts: How to restrict user to change (move) only one point in QGIS with Oracle Spatial database connected

texts:
I'm working with an Oracle Spatial View (Oracle Database 11g Enterprise, 64bit) in QGIS (2.18.0).

There was the bad case that I (and other user) was selecting a lot of points and accidentally moved them and saved the changes to the layer.

Now, I want to avoid this, with a database (trigger, constraint) or application (QGIS, PyQGIS) solution.

The user has to be allowed to move only one geometry feature at once. It should not be allowed to select two feature and move them.

How to distinguish between: 2 objects were changed at one time (not wanted) and objects were changed at different times which is ok.

My idea is to ask the user if the changes to the point geometry has to be saved after moving it. "Save Layer Edits" in QGIS has to disabled in some way and only the dialog saves the changes.

Or you restrict the user to select only one feature in this layer.

UPDATE - PyQGIS

The code checks how much features are selected, while the layer is in edit mode.

layer=None
for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "POINTS":
        layer = lyr

if layer.isEditable():
    selected_features = layer.selectedFeatures()

    selected_features_list = []
    for i in selected_features:
        selected_features_list.append(i["POINT_ID"])

    if len(selected_features_list) > 1:
        QMessageBox.warning(iface.mainWindow(),'Warning', 'You have selected more than one point. You can only move one point at once.')
        layer.removeSelection()
else:
    pass

This check has to be done between clicking/holding the select tool and releasing.

How can this be programmatically done?


texts: Dissolve Polygons using R / join specific polygons

texts:
I have data in shapefile format :: UK_boundaries
 I am new to R and Spatial data as well,  I am trying to join some areas to be as one area no line(s) between them. Based on my search I understand that  unionSpatialPolygons function dissolves polygons but when  I tried to use it  I got error:
Error: TopologyException: Input geom 0 is invalid: Self-intersection at or near point 162544.73551388 516922.55880385998 at 162544.73551388 516922.55880385998

libs <- c("rgdal", "maptools", "gridExtra", "rgeos", "raster")
lapply(libs, require, character.only = TRUE)

    LA <- readOGR(dsn = ".", layer = "infuse_dist_lyr_2011_clipped")
    IDs <- LA$geo_code

IDs[IDs %in% c("95DD","95II","95QQ","95EE","95RR","95UU","95TT","95CC","95BB")] <- c("95AA","95GG","95HH","95JJ","95JJ","95JJ","95KK","95OO","95XX")

LA_new <- unionSpatialPolygons(LA, IDs)

In summary I want areas of the same colour shown here to be joined.

Any help please?


texts: List attributes in field using ModelBuilder/ArcPy?

texts:
I'm using ArcGis 10.2.2 and I'm creating a model in ModelBuilder.

I want to know how I can list the attributes of a field and then make a selection for each unique attribute using Python.

I have a shapefile with the field "Name". Example:

name |
-----------
BB
BB
AB
AB
C
C

I want to list BB, AB, C, etc.

I want to construct a list and then use that list to make a selection for each of the unique values.  For example, "Name" = 'BB'.


texts: “What is the difference between ‘Select by location using a buffer or intersection’ and the ‘Buffer’ or ‘Intersect’ Analysis tools?


texts: Merged Layers turned into square shape cause by single layers in QGIS?

texts:
What I did is merged the 30 layers of the map (Isabela Map Province - Philippines) and below is what happened.

How do i remove the inside map of the border?

I'm new to QGIS and I'm using QGIS 2.4 Chugiak.


texts: Install GRASS GIS on macOS High Sierra 10.13 - error

texts:
I am trying to install for the first time GRASS GIS and after followed all the instructions from here https://grass.osgeo.org/download/software/mac-osx/ I tried to open the GRASS application and I got this:

wireless-staff-pt3-114-189:~ auuua$ '/Applications/GRASS-7.0.app/Contents/MacOS/grass.sh'; exit
Rebuilding Addon HTML manual pages index...
Rebuilding Addon menu...
Python 2.7.10 found.
Starting GRASS GIS...
dyld: Library not loaded: /usr/local/lib/libintl.8.dylib
  Referenced from: /Applications/GRASS-7.0.app/Contents/MacOS/lib/libgrass_gis.7.0.5.dylib
  Reason: image not found
Traceback (most recent call last):
  File "/Applications/GRASS-7.0.app/Contents/MacOS/gui/wxpython/gis_set.py", line 31, in <module>
    from core import globalvar
  File "/Applications/GRASS-7.0.app/Contents/MacOS/gui/wxpython/core/globalvar.py", line 29, in <module>
    from core.debug import Debug
  File "/Applications/GRASS-7.0.app/Contents/MacOS/gui/wxpython/core/debug.py", line 72, in <module>
    Debug = DebugMsg()
  File "/Applications/GRASS-7.0.app/Contents/MacOS/gui/wxpython/core/debug.py", line 37, in __init__
    self.SetLevel()
  File "/Applications/GRASS-7.0.app/Contents/MacOS/gui/wxpython/core/debug.py", line 43, in SetLevel
    self.debuglevel = int(grass.gisenv().get('WX_DEBUG', 0))
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 935, in gisenv
    s = read_command("g.gisenv", flags='n')
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 443, in read_command
    process = pipe_command(*args, **kwargs)
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 418, in pipe_command
    return start_command(*args, **kwargs)
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 356, in start_command
    if debug_level() > 0:
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 1519, in debug_level
    _debug_level = int(gisenv().get('DEBUG', 0))
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 935, in gisenv
    s = read_command("g.gisenv", flags='n')
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 446, in read_command
    return handle_errors(returncode, stdout, args, kwargs)
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 313, in handle_errors
    returncode=returncode)
grass.exceptions.CalledModuleError: Module run None ['g.gisenv', '-n'] ended with error
Process ended with non-zero return code -6. See errors in the (error) output.
Error in GUI startup. If necessary, please report this error to the GRASS developers.
Switching to text mode now.

Hit RETURN to continue...

Traceback (most recent call last):
  File "/Applications/GRASS-7.0.app/Contents/MacOS/grass70", line 1461, in <module>
    set_data()
  File "/Applications/GRASS-7.0.app/Contents/MacOS/grass70", line 682, in set_data
    gui_startup(grass_gui == 'gtext')
  File "/Applications/GRASS-7.0.app/Contents/MacOS/grass70", line 704, in gui_startup
    os.execlp(cmd_name, "-text")
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/os.py", line 329, in execlp
    execvp(file, args)
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/os.py", line 346, in execvp
    _execvpe(file, args)
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/os.py", line 382, in _execvpe
    func(fullname, *argrest)
OSError: [Errno 2] No such file or directory
logout
Saving session...
...copying shared history...
...saving history...truncating history files...
...completed.

[Process completed]

Any help?


texts: Combine multiple cached/tiled layers on a client-side openlayers application?

texts:
I have Geoserver with integrated geowebcache installed.

I have 3 layers: water, sewer, landbase, and each are cached separately.

Is it possible to toggle or combine all 3 cached layers into the openlayers application?
Is it efficient to combine all three separately like this? Would performance increase?


texts: Iterating through point shapefiles, converting them to polyline shapefiles with unique names using ArcPy

texts:
for fc in arcpy.ListFeatureClasses():
    outfc = fc + "_l"
    arcpy.PointsToLine_management(fc, outfc)

it comes up with the error "Output Feature Class: Same as input". Looks like it didn't convert at all, and added an "l" in the case of one file
Bacon.shpl
instead of what I want (Bacon_l.shp)

I am really bad with for loops and iterators, and I'd do this in ArcMap if I didn't have to manually enter a billion files using the batch feature.


texts: OpenLayers pass character parameters problem

texts:
I want send character parameter ('GREEN') to GeoServer as following:
var untiled = new ol.layer.Image({
    source: new ol.source.ImageWMS({
      ratio: 1,
      url: 'http://localhost:8080/geoserver/newyork/wms',
      params: {
            "LAYERS": 'newyork:viewStation',
            "viewparams": 'col:'+'GREEN'
      }
    })
  });

viewStation is a SQL query:
SELECT * FROM station WHERE color=%col%

but the code doesn't work, I also tried 'col:GREEN'
What can I change?


texts: Auto hide Attribute Table Widget in ArcGIS Viewers for Flex?

texts:
I'm trying to make a web application by using ArcGIS Viewers for Flex. I want that my Attribute Table can auto show up and hide by a button like happens on this site but I don't know how to do it.

Any suggestions?


texts: How to output geoJSON file from mysql spatial table?

texts:
I have a mysql spatial database ('mydb1') of one table ('pk'), i mean by spatial the fact that i obtained the database from a point shapefile using ogr2ogr tool. one row of this 'pk' table, is like:

 OGR_FID   SHAPE       CODIF    Position_X   Position_Y    Nom   Status
 1     [GEOMETRY-25o]  PAC182854  398235.38    414569.24    G-31   Vert

I need to output the geoJSON file from this spatial mysql table. I downloaded geoPHP library and also MySQL to GeoJSON script which i adapted to my database settings, like below:

<?php
/**
* Title:   MySQL to GeoJSON (Requires https://github.com/phayes/geoPHP)
* Notes:   Query a MySQL table or view and return the results in GeoJSON format, suitable for use in OpenLayers, Leaflet, etc.
* Author:  Bryan R. McBride, GISP
* Contact: bryanmcbride.com
* GitHub:  https://github.com/bmcbride/PHP-Database-GeoJSON
*/

# Include required geoPHP library and define wkb_to_json function
include_once('geoPHP/geoPHP.inc');
function wkb_to_json($wkb) {
    $geom = geoPHP::load($wkb,'wkb');
    return $geom->out('json');
}

# Connect to MySQL database
$conn = new PDO('mysql:host=localhost;dbname=mydb1','root','admin');

# Build SQL SELECT statement and return the geometry as a WKB element
$sql = 'SELECT *, AsWKB(SHAPE) AS wkb FROM pk';

# Try query or error
$rs = $conn->query($sql);
if (!$rs) {
    echo 'An SQL error occured.\n';
    exit;
}

# Build GeoJSON feature collection array
$geojson = array(
   'type'      => 'FeatureCollection',
   'features'  => array()
);

# Loop through rows to build feature arrays
while ($row = $rs->fetch(PDO::FETCH_ASSOC)) {
    $properties = $row;
    # Remove wkb and geometry fields from properties
    unset($properties['wkb']);
    unset($properties['SHAPE']);
    $feature = array(
         'type' => 'Feature',
         'geometry' => json_decode(wkb_to_json($row['wkb'])),
         'properties' => $properties
    );
    # Add feature arrays to feature collection array
    array_push($geojson['features'], $feature);
}
header('Content-type: application/json');
//write to json file
$fp = fopen('pk_geojson.json', 'w');
fwrite($fp, json_encode($geojson));
fclose($fp);
$conn = NULL;
?>

But when i execute the code on the browser and open the output geojson file, i don't find the longitude and latitude values of each feature. i get instead, X and Y values another time, like so:

{"type":"FeatureCollection","features":[{"type":"Feature","geometry":{"type":"Point","coordinates":[398235.3882,414569.2462]},"properties":{"OGR_FID":"1","codif":"PAC182854","position_x":"398235.388200","position_y":"414569.246200","nom":"G-31","point_kilo":"182000","status":"Vert"}}]}

in the geometry property like you can see, i get {"type":"Point","coordinates":[398235.3882,414569.2462]} so another time x, y values that i have in the property "properties" ("position_x", "position_y"), while i need to have instead longitude and latitude values, in order to add this geojson file to my map on openlayers as vector layer. I notice that in the starting, i had a x,Y excel file that i converted with arcmap to the shapefile from which i built mysql database using ogr2ogr tool. So, is it a problem of the way i made the shapefile or what's wrong exactly i have to fix to get a proper geojson file?


texts:
WMS standard allows specifying a list of styles for each WMS layer. The users can then use the style name in the URL query parameter "STYLES". If the styles are not specified, "default" is assumed.

But, how do I select such a style in QGIS? I was able to add a WMS server, and add and display its layers. However, although I can see all the styles listed in the "Layer properties", I haven't find any option how to select one.

For example, if some WMS service has a layer named "Vegetation", and for this layer the service defines styles "Grayscale" and "Green", how do I select "Green" instead of the default one?


texts: How to polygonize using GeoDjango and PostGIS?

texts:
The code of GeoDjango-1.9.2 contains reference to the function ST_Polygonize at:
https://github.com/django/django/blob/master/django/contrib/gis/db/backends/postgis/operations.py#L132

But is it seems like not documented anywhere.
Is the function available for use somehow?


texts: Network Analyst Add City Boundary as Polygon Restriction

texts:
I have a network I am trying to create but am finding the service area polygon I create are "leaking" outside of the city boundary.

I was curious if there was a way to add the polygon as a restriction so that these polygon are "snapped" to the city polygon?

I have tried to add it as a restriction polygon but I come up with errors for the multiple locations I am attempting to create a service area for. I also converted my polygon to a line and utilized the snap feature in Layer Properties, and that didn't work either.