texts:
I am trying to write a Python command that includes a logical query. I have series of tables and layers and would like to run a command on them when two shared variables have the same values. I created a list of those tables and layers:
pattern = '*Assess*.dbf'
files = []
tables = []
layernames = []
## create the loop element
for dirpath, dirnames, filenames in arcpy.da.Walk(inputFolder):
    # append shapefile path to list
    for filename in fnmatch.filter(filenames, pattern):
        files.append(os.path.join(dirpath, filename))

tables = arcpy.ListTables()

for file in files:
    names.append("\"TOWN_ID\" = " + str(int(os.path.basename(file).rstrip(os.path.splitext(file)[1])[1:4])))
    layernames.append(str(int(os.path.basename(file).rstrip(os.path.splitext(file)[1])[1:4])))

The code is as followed:
for layername, table in zip(layernames, tables):
    #print layername
    #print table
    #CF_QueryTool = "C:\\ArcGIS_Temp\\parcel_mosaic\\land_parcels_mosaic.gdb\\cf_QueryTool"
    tablelist = [table, "TaxParcels_Layer.{}".format(layername)]
    whereclause = "TaxParcels_Layer.{}.LOC_ID = table.LOC_ID".format(layername)
    print tablelist
    print whereclause
    arcpy.MakeQueryTable_management(tablelist, "QueryTable.{}".format(layername), "ADD_VIRTUAL_KEY_FIELD", "", "" ,whereclause)
    print "QueryTable.{} created".format(layername)

My where clause translates to "TaxParcels_Layer.1.LOC_ID = table.LOC_ID".
The problem is, apparently there is something wrong with my WHERE-clause which made me wonder if I can actually use something like table.locid using a list of tables. Let me know if I'm making a mistake here.


texts: Finding outliers with wrong Z values in LIDAR data

texts:
I have more than 2000 XYZ (only ground points) .txt files and I have to do Quality Control (QC).

Let's say I have file, which I converted to multipoint (or even singlepart feature). I have mean Z value for a file, but if most of the points are 20-40 metres, and two, three points per file are 70-80 metres (wrong Z values). How to catch these points and remove? On another hand, it could be that in another file 70-80 metres is common Z value and some points wrong Z value 120 metres. In another file 120 metres is common and correct.

How to find and remove points with wrong Z values in LIDAR data?

Standard deviation? Raster cells?


texts: TileMill Missing Streets HERE street data

texts:
I am creating tiles with TileMill. Before adding layers to TileMill I'm checking the geometries for errors and I'm not seeing any

However, once I add the street layer to TileMill there are missing streets

I thought maybe this was isolated to a rendering issue in TileMill and would disappear after I generated tiles, but the problem persisted afterwards when I loaded it into GeoServer.

This is the styling I'm using in TileMill

#MtViewStreets3CLs[zoom>12]{
line-width:1;
line-color:#b6b6b6;
}

#MtViewStreets3CLs[zoom>14]{
  line-width:1;
  line-color:#b6b6b6;
  text-name:[ST_NAME];
  text-face-name:"Arial Regular";
  text-fill:#838383;
  text-size: 11;
  text-placement: line;
  text-min-path-length:100;
  text-avoid-edges:true;
  text-min-distance:50;
  text-dy: 6;
  text-max-char-angle-delta: 15;
}

I checked this post but it's not the same issue I'm having:

Tilemill didn't render streets correctly, some parts are missing


texts:
How do I add a location attribute to a point file, where the location is an ID from an underlying shapefile, using ArcMap 10.1?

Two layers, I have polygons with a location ID, 1,2,3,4,5 etc. I have point layer which I would like to add a new column to, and update that column with the ID from the polygon file


texts: ST_ConcaveHull interaction with ST_Boundary


texts: Crop a raster with another raster in R / GIS

texts:
I have raster files with the same resolution (30x30m), but different dimensions, and I would like to have all rasters with the same dimension and resolution.

So, here is the first raster:

(A)

and here is the second:

(B)

Together, they look like this:

I want them to have the same shape. Using R, I managed to crop A according to B, by simply using the raster::crop() function, and created a new raster (let's call it A1). However, it doesn't work if I try to crop B according to A or A1, but I don't get any kind of error. I am even able to plot it, and it shows the raster with the same rectangular shape, instead of considering the limits of A/A1...

I have also tried with QGIS and I have exactly the same problem: I can crop A according to B but not the opposite.

Does anyone know how can I crop this raster?


texts: PostGIS query with input from point layer

texts:
I have a table of rivers in PostGIS, each river linking to its outflow.
I also have a QGIS layer of some selected sources, and would like to select the entire river downwards from each source (i.e. find the closest segment to each point, then recurse through the outflow links until hitting the ocean).

Is there a way to do this interactively? I'd like to just update the point layer and have the query use these points. I don't think QGIS offers any help with that (parametrized queries?), but I could store the point layer as a shapefile, load that as a temporary table and use that in my query.

is there a way to load a shapefile in a subquery? I don't want to use shp2pgsql manually, but load the current version of the file each time the query is run.
can I maybe somehow trigger actions in QGIS when a layer has been edited? i.e. update the river layer when the point layer has been saved with the new query output, or, alternatively, run a script that loads the point data, runs the query, dumps its data into a shapefile, loads that in QGIS and refreshes the view?


texts:
I can't find a way to clear away all imported QGIS styles, even completely removing all QGIS files, after new install - old styles are still there. Where to look?

Removing one by one is the only option from within the application, but there are hundreds that I have imported.

I found the "first run state" question, but that answer works for Windows, but how about MacOS?


texts: Inserting geography data types with GeoPandas

texts:
I would like to insert a geodataframe containing point geometries in SRID 4326 to an existing table in PostGIS+PostgreSQL.  The issue is that I would like to retain the existing table structure, which has GEOGRAPHY data type column, but inserts fail.  Using a sqlalchemy connection engine, I try the following:
dtypes(gdf)
ts              datetime64[ns]
geog                  geometry
id                      int64

dtype: object

The geodataframe has a geometry column named geog
import geopandas as gpd
import psycopg2
import numpy as np
import json
from collections import defaultdict
import sqlalchemy
from geopandas import GeoDataFrame
from shapely.geometry import Point
from geoalchemy2 import Geometry, WKTElement, Geography

with engine.connect() as conn:
    gdf.to_postgis('my_data', con=conn, if_exists='append', index=False, dtype={'geog': Geography('POINT', srid=4326)})
    conn.close()

but this returns an error:

 (psycopg2.errors.RaiseException) find_srid() - could not find the corresponding SRID - is the geometry registered in the GEOMETRY_COLUMNS table?  Is there an uppercase/lowercase mismatch?
CONTEXT:  PL/pgSQL function find_srid(character varying,character varying,character varying) line 17 at RAISE

[SQL: SELECT Find_SRID('public', ' my_data', 'geog');]

The column names and types match exactly and in PostgreSQL, the geog column is of type geography:
CREATE TABLE public.my_data (

    id   SERIAL PRIMARY KEY,
    geog GEOGRAPHY(Point, 4326) NOT NULL,
    ts TIMESTAMP NOT NULL

)

What is the correct way to insert geodataframes with GEOGRAPHY data type into PostgreSQL?


texts: Esri Leaflet refresh() method not working

texts:
I'm looking to refresh a featureLayer in Esri Leaflet with the refresh() method (API reference: https://esri.github.io/esri-leaflet/api-reference/layers/feature-layer.html).

I set up a button to refresh the featureLayer:

const testLayer = L.esri.featureLayer({url:''}).addTo(map);

document.getElementById("test").onclick = function(){
testLayer.refresh();
});

However, I would like to set it up with a setInterval method. So when people are making changes and editing data it will update the webmap dashboard.

setInterval(function() {
    testLayer.refresh();
}, 3000);

I need the feature layer to refresh so when there's changes to the data that webmap will reflect that without reloading the page. The ArcGIS API for Javascript has a similar and working method documented here support.esri.com/en/technical-article/000013384 -- however I'm using the esri Leaflet plugin. It seems to pretty simple to me but I'll set it up and make a change to the data and the map will not update the featureLayer unless I press refresh on the page.


texts:
I have two tablets with QField on it, on which I load the same QGIS Project. The Fieldworkes go record points in the field. After that I want to synchronize both tablets back to my QGIS Project via QSync. After the first synchronization I've got the message: "CRITICAL Data from this file are already synchronized with the original project".
Both the tablets have the same data on it but the fieldworkers record the points in different areas.

Does anyone have an idea how I can syncronize two tablets?


texts: Can anyone see what I did wrong with this Make Feature layer function? (Arcpy)

texts:
So I am just testing this function by itself before I put it in my loop. So I run the test (Python 2.7) and I get a result saying that it was create, but I do not see it in folder I assigned it to. I even tried it on the arcpy box in arcgis 10.2 and I still can't find the layer. I double checked the path I directed it too but nothing.

    import arcpy
    import os

    arcpy.env.workspace = 'C://Users/my/workspace/'
    arcpy.env.overwriteOutput = True

    inPut = 'C://Users/my/workspace/'

    allStorms = 'allstorms.shp'

    acrpy.MakeFeatureLayer_management (allstorms, inPut+'group1')


texts:
My question is simple and hopefully can be answered quickly. I have a code in R that is for calculating home ranges. Everything works except one line of the code, which I've included here:

locoh <- lapply(dat, function(x) rhrLoCoH(x[ c("LONGITUDE", "LATITUDE")]))

The error I get is as follows:

Error in `[.data.frame`(x, c("LONGITUDE", "LATITUDE")) : undefined columns selected

My guess is that I have messed up somewhere with my commas, but I'm not sure where or how to fix it.


texts:
I am trying to import a shapefile using db manager to postgresql. My shapefile has OBJECTID (PRIMARY KEY), ID (STRING FIELD eg: 2543-96-7948 )and several other attributes of type string.
When importing I get the following error.

When I import using postgis shapefile import/export manager, I got no errors and the shapefile uploaded successfully.
I think qgis is assuming my primary key to be ID field when it is actually OBJECTID. How to change this?
I tried using the primary key field in input table and set it to id but got same error. Does it create new field or is it selecting a column from the attribute table?
looked up for answers here but was not enough to resolve my issue.


texts:
I am using QGIS 2.16.
I have a database of 100 entries which has been categorised using a column of data that has 20 categories which essentially breaks up my map into 20 regions.
How do I create an automatic  boundary around the 20 sub groups?

I went down the voronoi and dissolve route, which I think works but I have some queries to work out kinks. I'm unsure if this should be a new question entirely.
In points 2,3,4 there are unwanted polygons.
How do I remove these?
Whilst in point 1, the symbol within the triangle should black as it is a different category.
How do i fix this?
Is it my 20% buffer that is causing the issue?
Finally, how do I enlarge the 'frame' around the map of Scotland as it cuts through the Isle of Orkney (top of the map)?


texts: Tips for learning to work with geoserver in Javascript

texts:
After succesfully being able to work with the ArcGIS Javascript API in a webpage I want to learn to work with Geoserver. I have downloaded Geoserver and uploaded one shapefile. Now I am trying to call this shapefile from a Javascript/HTML script. For the ArcGIS Javascript API there is a lot of documentation. For geoserver a lot less is to be found. All I found is this: http://docs.geoserver.org/2.6.x/en/user/community/scripting/js/api/index.html#geoserver-javascript-api-documentation

I am not succeeding with just that. Could any of you give me a source or advise to learn more about working with geoserver in a webpage?


texts: Copying the numeric portion of an address to a new field

texts:
I'm attempting to copy just the numeric portion of a street address into a new table using the Field Calculator in QGIS. The expression I'm trying to use is:

left(  "SITE_ADDR", strpos(  "SITE_ADDR", " " ) -1)

strpos is giving an invalid expression error- column " " not found.

Is there another way to format this expression to actually work?


texts: Test for (curly) brackets in expressions in QGIS 2.14.14 with strpos() weird results

texts:
When trying to find a strpos() of a substring containing curly brackets (or brackets anyway as far as I can see) I get 0 as a result.

Test for substring without {} works fine:

Test for substring with {} weird result (should be 6):

Finding the brackets seems no problem, because replace() works as expected:

Aren't brackets treated as normal characters in strpos(), am I missing something? Bug? Feature?


texts: Running viewshed analysis in ArcGIS Online

texts: Can I run a viewshed analysis in ArcGIS Online, or would I need an ArcGIS Desktop / Spatial Analyst licence?


texts: Interpreting the output tables from point distance measures from two coordinates in ArcGIS

texts:
I have recently used ArcGIS and prior to this I have no experience with it. I wanted to calculate the distances between Care Homes in England to NHS hospitals. Basically I wanted to see how many Care Homes were there within a radius of 10 km to each hospital. I used the point distance method and in the options menu in the search radius option I specified 10 kilometers.

This is the resultant table that I get.

My query is the distance column are expressed in all decimal representations and therefore I don't think they are in km. For example, the distance between INPUT_FID 0 and NEAR_FID 20168 is 0.08618273742. In terms of kilometers, the distance between their coordinates however is 7.880

Do you have any explanation as to why this discrepancy arises?


texts: Setting an identical map extent for all systems resolutions with OpenLayers

texts:
I want to see my map extent as identical in all resolutions (17 inches, 20 inches etc).

If I set centre and zoom together, the map extent is varying with respect to different machine resolutions.

view: new olView({
            center: [312234.8270497762,2636116.2213047906],
           // extent: mapExtent,
            projection: this.projection,
            rotation: 0,
            zoom: 8
          })

How can I fix it?


texts: ArcGIS Server won’t release SDE table locks

texts: My ArcGIS Server services are locking tables in SDE. This is a problem if I need to edit or overwrite them. Unfortunately the locks don't go away upon stopping the services. The only solution I have found is to restart the Linux arcsde service, which has many implications. Is there a better way to lose these phantom locks?


texts:
Is it possible to generate a kml file of searched items in Google Maps. For example:

Search for highschools in Virginia
Return the results as a kml file with the schoolname address and location


texts: Displaying attribute field value in TextElement gives AttributeError: 'list' object has no attribute 'text'?

texts:
I am new in Python scripting. I have a geodatabase with layer Test and this layer has field Name and Text element on layout named TEname. I am trying to show selected feature attribute Name value in Text element, but I could not do it.  I can print the value for selected feature, but can't show it in Text element.

Here is the code I use:

import arcpy

mxd = arcpy.mapping.MapDocument("current")
TEname = arcpy.mapping.ListLayoutElements(mxd,"TEXT_ELEMENT","Name")[0]

for row in arcpy.da.SearchCursor("Test","name")
     print row

...
(u'Land plot 1',)

for row in arcpy.da.SearchCursor("Test","name"):
    TEname.text = row

...

      TEname.text = row

Runtime error
Traceback (most recent call last):
  File "", line 1, in
  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\arcobjects_base.py", line 89, in _set
    return setattr(self._arc_object, attr_name, cval(val))
RuntimeError: TextElementObject: Error in setting text


texts: Specify Region for Export Table in Google Earth Engine?

texts:
I have a dataset of polygons for the entire planet that I'd like to export into a KML/KMZ file. However, I only need the data for South America. Is is possible to specify which regions to export?
This is the code I currently have:
var polygons = ee.FeatureCollection('WCMC/WDPA/current/polygons');

Export.table.toDrive({
  collection: polygons,
  folder: 'GEE_exports',
  fileFormat: 'KMZ'
});


texts: How to get the name of the layer on a map composer with a single layer?

texts:
I tried

layer_property( 'layer', 'title')

I would like to return the current layer from the map instead of 'layer'.
(or list of the layers on the map?)

Doesn't seems to works with

 layer_id

or

@layer_name

It's for a model composer, that's why i want it dynamics.


texts: Accessing ArcView 3.X help after Windows 7 install?

texts:
How do I access ArcView 3.X help after Windows 7 install?

ArcView works but I can't access help.


texts: looking for a GPS device which can transmit the location directly

texts:
I am looking for a GPS device or receiver that can transmit its own location to a specified configured server directly.

there will be no medium like mobile phone or computer in between GPS and Server which receives the data.

a dongle sort of thing which has a sim card in it and with the use of that it can transfer data directly to server


texts: sorry if this is a simple question, but based on my search terms I keep getting PDF papers about DEMs. Someone has asked me to generate fill quantities off of their lidar data which they have presented to me as a PDF. I'm under the assumption that a  PDF is treated as a flat image and can not produce elevation information to a GIS program - but am I wrong? I converted the PDF to a .tiff composed of 4 bands 0-255 which screams RGB to me....


texts: Vector Geoprocessing Tools Missing from QGIS 2.18.2 for Mac OSX 10.11.6

texts:
I tried the steps from "Vector Tools Missing in QGIS 1.8 Ubuntu".
 Most of the tools under "Vector" are missing from my recent installation of QGIS 2.18.2, so I am unable to clip, buffer, make Voronoi polygons, etc. Python appears to be running correctly and I have checked to ensure my plugins are updated and enabled. For some reason "ftools" plugin is not availabe when I search for it in Plugin Manager under "All".

EDIT: I did a clean reinstall of QGIS and this fixed the problem. I think there was an old version of GDAL causing the problem.


texts:
I would like to render contours using Mapnik. The problem is that some contours are very long which makes rendering too slow and so I decided to spling long contours to multiple smaller segments. Input table (generated with gdal_contour):

gis=> \d contour;
                                            Table "public.contour"
    Column    |            Type             | Collation | Nullable |                 Default
--------------+-----------------------------+-----------+----------+------------------------------------------
 ogc_fid      | integer                     |           | not null | nextval('contour_ogc_fid_seq'::regclass)
 id           | numeric(8,0)                |           |          |
 height       | numeric(12,3)               |           |          |
 wkb_geometry | geometry(LineString,900914) |           |          |
Indexes:
    "contour_pkey" PRIMARY KEY, btree (ogc_fid)
    "contour_wkb_geometry_geom_idx" gist (wkb_geometry)

I am no GIS or PostGIS proffesional but I tried my best and came up with following query:

create table contour_split as select
  ogc_fid,
  height,
  (select st_makeline((x).geom) from (select st_dumppoints(wkb_geometry) as x offset round(start) limit round(start + seg_len) - round(start)) as sub3) as geom
from
  (
    select
      ogc_fid,
      wkb_geometry,
      height,
      n_pts / ceil(n_pts / 1000.0) as seg_len,
      generate_series(
        0,
        n_pts - n_pts / ceil(n_pts / 1000.0),
        n_pts / ceil(n_pts / 1000.0)
      ) as start
    from
      (select *, st_numpoints(wkb_geometry) as n_pts from contour) as sub1
  ) as sub2;

The problem is that this query is very slow. There are over 1.6M rows in the contours table and the query is running for several hours already.

I've also tried to do it with intermediate tables but the speed (of the last query) is still slow:

alter table contour add column np int;

update contour set np = st_numpoints(wkb_geometry);

create table contour_tmp as select
    ogc_fid,
      generate_series(
        0,
        np - np / ceil(np / 1000.0),
        np / ceil(np / 1000.0)
      ) as start
    from contour;

create table contour_split as
select
  ogc_fid
  height,
  (
    select
      st_makeline((x).geom)
    from
      (
        select
          st_dumppoints(wkb_geometry) as x offset round(start)
        limit
          round(start + np / ceil(np / 1000.0)) - round(start)
      ) as ddd
  )
from
  contour
  join contour_tmp using (ogc_fid);

Can somebody see there something "obvious" what could improve the speed of the query? Is there a better approach to achieve even splitting linestrings to not to be longer than 1000 points? ST_Subdivide returns uneven results - some segments are very short. Or to make Mapnik fast even with long contour linestrings?

Version information: POSTGIS="2.5.1 r17027" [EXTENSION] PGSQL="110" GEOS="3.7.1-CAPI-1.11.1 27a5e771" PROJ="Rel. 5.2.0, September 15th, 2018" GDAL="GDAL 2.4.0, released 2018/12/14" LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.3.1" RASTER


texts: I don't want to use the time_of_day parameter in MakeServiceAreaAnalysisLayer, what do I put to miss it?

texts:
I'm modelling accessibility, so am using network analyst arcpy.na.MakeServiceAreaAnalysisLayer().

However, I don't want to enter anything for the time_of_day parameter. How do I write my Python script so that I can write parameters after this, but not include time_of_day or time_zone?

I get ERROR 000628: Cannot set input into parameter time_of_day.

I'm a newbie to scripting.


texts: Understanding precision of ArcGIS geodatabases?

texts: I would like to know how precision is handled in ArcGIS software. Is it using a real Computer Algebra System engine, an arbitrary precision library or do others tricks to handle precision ?


texts: I'm wondering if it is possible to adjust the cell height of the attribute table (QGIS 2.14.3). I would like to reduce the cell height, to avoid the unused space in a cell.


texts:
I need to create a dialog that the user can select a proper Workspace.

Is there a similar to IAddDataDialog, but for workspaces?


texts: Magnetic declination - which magnetic model to choose?

texts: which model to determine the magnetic declination is more accurate WMM (2014-2019) or IGRF (1590-2019)? I just wanted to determine the magnetic declination for calibration TruPulse Laser Rangefinder. I did not know which model would give me greater precision measurements (Azimuth+Distance)+(GPS+GIS)


texts: Strange results with Voronoi polygons in QGIS


texts: Adding date field to attribute table using QGIS processing modeler?

texts:
I want to add a date field to my output layers attribute table, using processing modeler. I want the Date field to be a parameter in my model, because I want to manually input the  Date everytime I run my model. Is there a way to do this, so the Date appears as a parameter before running the model (as shown in the picture) and the date entered then appears in the attribute table under a field titled Date. I want to give each layer I create a Date in the attribute table without having to edit the model each time.

I am new to QGIS.


texts: How to find name of tables for which you are the owner in pgAdmin 4?


texts:
I'm trying to make a sum of n rasters and create a new one using QgsRasterClaculator. The problem is, every time the code reaches processCalculation, it makes QGIS close and create a dump file. I don't know if I'm doing it wrong, but I'd like to ask you for some help in how I can do it.

That's how I'm trying to do:

entries = []
rasters = os.listdir("C://test/Rasters/")
for i in range(0, len(rasters)):
        raster = "C://test/Rasters/"+rasters[i]

        resp = QgsRasterLayer(raster)
        rast1 = QgsRasterCalculatorEntry()
        rast1.raster = resp
        rast1.ref = "rast"+str(i+1)
        print(rast1.ref)
        rast1.bandNumber = 1
        entries.append(rast1)

        if(i == len(rasters)-1):
            var = var + 'rast'+str(i+1)
        else:
            var = var + 'rast'+str(i+1)+' + '

calc = QgsRasterCalculator( var, 'C://test/Rasters/raster_somado.tif', 'GTiff', resp.extent(), resp.width(), resp.height(), entries )
print(calc.processCalculation())

I think the error might be happening because I'm using a variable (var), that contains the rasters alias. The base code I found and that works is

calc = QgsRasterCalculator( 'rast1 + rast2', 'C://test/Rasters/raster_somado.tif', 'GTiff', resp.extent(), resp.width(), resp.height(), entries )
print(calc.processCalculation())

Can anybody help me? I already tried some things but it's still not working, QGIS keeps creating the dump and closes.


texts: Python Labeling using dictionary not printing all values

texts:
I am labeling manhole inverts using python and ArcMap 10.3. When I use a dictionary to control the label order some manholes don't print all of the inverts while other manholes do. I can post field information/attribute table info if needed. I suspect its an error on my part using an empty/none field in the dictionary and underlying hash table but I really don't know enough to say that with confidence.

How can I get all the inverts to display? - do a check on every argument for none or "" before inserting into the dictionary?

Working Label Expression:

I think here dict is set to a list with grouped values (key,value),(key,value),etc and everything is working correctly:

import json
def FindLabel ( [Elevation] , [Out_Elev] , [Out_Loc] , [In1_Elev] , [In1_Loc] , [In2_Elev] , [In2_Loc] , [In3_Elev] , [In3_Loc]  ):
  dict = ([Out_Elev],[Out_Loc]),([In1_Elev],[In1_Loc]),([In2_Elev],[In2_Loc]),([In3_Elev],[In3_Loc])
  x = json.dumps(dict, sort_keys=True)
  return x

Not-Working Label Expressions:

The next two pictures show where I try to use a dictionary. This first one is changing dict from a list to a dictionary using dict = OrderedDict(dict) and the second one is setting straight to a dictionary using dict = {key:value,key:value,etc}. I don't understand why all 3 inverts display for the manhole at the top of the picture but not for the one on bottom left.

Using dict = OrderedDict(dict)

from collections import OrderedDict
import json
def FindLabel ( [Elevation] , [Out_Elev] , [Out_Loc] , [In1_Elev] , [In1_Loc] , [In2_Elev] , [In2_Loc] , [In3_Elev] , [In3_Loc]  ):
  dict = ([Out_Elev],[Out_Loc]),([In1_Elev],[In1_Loc]),([In2_Elev],[In2_Loc]),([In3_Elev],[In3_Loc])
  dict = OrderedDict(dict)
  x = json.dumps(dict, sort_keys=True)
  return x

Another code variant where inverts are not showing using dict = {key:value,key:value,etc}

import json
def FindLabel ( [Elevation] , [Out_Elev] , [Out_Loc] , [In1_Elev] , [In1_Loc] , [In2_Elev] , [In2_Loc] , [In3_Elev] , [In3_Loc]  ):
  dict = {[Out_Elev]:[Out_Loc],[In1_Elev]:[In1_Loc],[In2_Elev]:[In2_Loc],[In3_Elev]:[In3_Loc],}
  x = json.dumps(dict, sort_keys=True)
  return x


texts: How to find the nearest polygon from a given point?

texts: Basically I have set of points [lat-long pairs] which define different polygons on a map. Now I have a coordinate [lat-long] from which I want to find the nearest polygon from that point.


texts: Store Geometry in a database with the ESRI Silverlight API

texts:
Is it possible to store geometry objects in a database using the ESRI Silverlight API?  For this particular case, the client may have SQL Server 2005 installed instead of 2008 and therefore I can not use spatial types.

Any solution would also have to support clients who use ArcGIS Server 9.3, and therefore we would not be able to use the new editing capabilities in ArcGIS 10.


texts: R: Substracting rasters by their minimum values (loop over a list)

texts:
I'm trying to reclassify a bunch of rasters by substracting the pixels of each of them with their minimum value (eg: I want to substract all the pixels of a raster(A) by its minimum value, same thing for raster(B), (C) etc...).
It was not a problem to do such task in QGIS when I only had a few ones but now that I have more than 200 rasters I don't have other choice than automating the process.
Here is my script so far:
EDIT:
library(raster)

mypath <- "F:/Stage/Modele_SIG/Etape_2/test/"

files <- list.files(mypath, pattern="*tif$")

rasters <- lapply(paste0(mypath, files), raster)

for (i in 1:length(rasters)) {
  r <- rasters[[i]]
  r <- r - cellStats(r, "min")
  r <- writeRaster(r, files[i], overwrite=TRUE)
}

r

Unfortunately it seems that only one of my rasters is fetched in rasters[i] when the loop is run whereas I have all my rasters in rasters
Any ideas about what I'm doing wrong and how I could make my script works ? Thanks in advance !


texts: How to get attribute values from snapped feature in QField

texts:
Is there anyway to get attribute values automatically filled from snapped feature when digitizing a new line?
E.g I have a previous line digitized with attribute "Code" and its value "100". Then I would like to digitize a new line starting from snapped node of that previous line with same attribute values prefilled in list before saving the feature.
Just to clarify, QField's pcre-built option to keep last used value in list is not what I want.


texts:
I have a .mat file I got from someone, and it's supposed to represent global gridded data. I wanted to read it as a raster in R.

Importing was ok, I was able to import using the R.matlab package

> library(R.matlab)
> library(raster)
>
> data <- readMat("map.mat")
> str(data)
List of 1
 $ comData: num [1:360, 1:720] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
 - attr(*, "header")=List of 3
  ..$ description: chr "MATLAB 5.0 MAT-file, Platform: PCWIN64, Created on: Mon May 11 09:17:06 2015                                        "
  ..$ version    : chr "5"
  ..$ endian     : chr "little"

I would have expected at least three lists, with longitude, latitude, and the data. However, it appears to only have one list comData. I have never imported data from Matlab, I am fairly new with GIS, so probably I am missing something. How can I read this dataset as a raster?


texts: Calculating distances between points two feature datasets

texts:
I will do my best to explain my problem. I have two feature datasets. One set has about 100 points that are "city centers." The other dataset has points for each census tract in each of these 100 cities. THere is a row of data per census tract as well as an identifier number for each city (which is also in the city center feature dataset).

I need to measure the distance between each city center and each point/census tract in that city and then have these distances in a table with the census tracts numbers....I was hoping there was some way to use python (I'm a beginner at python!) to loop through each point and measure the distance to the city center.

ANY suggestions on how to best approach this would be great and the more detailed your response, the better! I've used point distance tool before when I only had 5 cities so I could just run each one manually to get the distances (each city was also a separate file)....


texts: How to remove extra field from LAS file with laspy

texts:
I have LAS files with a useless extra scalar field (named "tree") that I would like to remove.
I tried the following code to create a new LAS file, copy each field of the orginal LAS file into the new LAS file, except the unwanted field.
import laspy

LASInput = laspy.file.File('foo.las', mode = 'r')
LASOutput = laspy.file.File('bar.las', mode = "w", header = LASInput.header)

#Lets copy all the existing data from input file:
for dimension in LASInput.point_format:
    print(dimension.name)
    if not(dimension.name == 'tree'):
        LASOutput.writer.set_dimension(dimension.name, LASInput.reader.get_dimension(dimension.name))

LASOutput.close()
LASInput.close()

But the resulting LAS file still contains the field "tree". I guess the header also has to be modified so I tried this approach, creating a new header:
import laspy

LASInput = laspy.file.File('D:\\foo.las', mode = 'r')
LASOutput = laspy.file.File('D:\\bar.las', mode = "w", header = laspy.header.Header(1.3))

#Lets create the extra fields
NamesOutput = [i.name for i in LASOutput.point_format]
for dimension in LASInput.point_format:
    print(dimension.name)
    if not(dimension.name == 'tree'):
        if not(dimension.name in NamesOutput):
            LASOutput.define_new_dimension(name = dimension.name, data_type = 5, description ='blabla')

#Lets copy all the existing data from the imput file
for dimension in LASInput.point_format:
    print(dimension.name)
    if not(dimension.name == 'tree'):
        LASOutput.writer.set_dimension(dimension.name, LASInput.reader.get_dimension(dimension.name))

LASOutput.close()
LASInput.close()

print('ok')

This works, but I find it rather inelegant, and also I did not find a way to get the fromat of the fields in the input header, in my code I used data_type=5, which is likely to be incorrect. Is there not a way to duplicate the input header and then to discard the unwanted field from the header?


texts:
When running cmake for CGAL 3.9, I get the following output:

root@postgresql local/CGAL-3.9# cmake -DWITH_CGAL_Qt3=OFF -DWITH_CGAL_Qt4=OFF -DWITH_CGAL_ImageIO=OFF .
-- The C compiler identification is GNU
-- The CXX compiler identification is GNU
-- Check for working C compiler: /usr/bin/gcc
-- Check for working C compiler: /usr/bin/gcc -- works
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Check for working CXX compiler: /usr/bin/c++
-- Check for working CXX compiler: /usr/bin/c++ -- works
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Build CGAL from release: CGAL-3.9
-- Packagenames: CGAL-3.9
-- CGAL_MAJOR_VERSION=3
-- CGAL_MINOR_VERSION=9
-- CGAL_BUGFIX_VERSION=
-- CGAL_SONAME_VERSION=8
-- CGAL_SOVERSION     =8.0.0
-- CGAL_REFERENCE_CACHE_DIR=
-- Building shared libraries
-- Targetting Unix Makefiles
-- Using /usr/bin/c++ compiler.
-- USING CMake version: 2.8.0
-- System: Linux
-- USING GCC_VERSION = '4.4.3'
-- Using gcc version 4 or later. Adding -frounding-math
-- Build type: Release
-- USING CXXFLAGS = ' -frounding-math -O3 -DNDEBUG'
-- USING LDFLAGS = ' '
-- Could NOT find GMP
-- Could NOT find MPFR
-- CGAL_Core needs GMP, cannot be configured.
CMake Error at cmake/modules/FindBoost.cmake:917 (message):
  Unable to find the requested Boost libraries.

  Boost version: 1.40.0

  Boost include path: /usr/include

  The following Boost libraries could not be found:

          boost_thread

  No Boost libraries were found.  You may need to set Boost_LIBRARYDIR to the
  directory containing Boost libraries or BOOST_ROOT to the location of
  Boost.
Call Stack (most recent call first):
  cmake/modules/CGAL_SetupBoost.cmake:12 (find_package)
  cmake/modules/CGAL_SetupDependencies.cmake:32 (include)
  CMakeLists.txt:422 (include)

-- Boost include:     /usr/include
-- Boost libraries:
-- Boost definitions:
-- USING BOOST_VERSION = '1.40.0'
...

I've tried various env variables.  I am wondering if I need to install / compile some other boost thing?  I have some libboost* in /usr/lib
Any suggestions on things to try would be appreciated.


texts: How to export CSV data of area from Google Earth

texts:
I want to export CSV data of some area in Google Earth.
I will then import this data into python and run a 3d map.
I will be able to rotate and scale it then.

I can do the last 2 steps, but how can I do the 1st step?

I want to export it like that with that code:

import plotly.graph_objects as go

import pandas as pd
# here read data:
z_data = pd.read_csv('https://raw.githubusercontent.com/plotly/datasets/master/api_docs/mt_bruno_elevation.csv')
# plot it:
fig = go.Figure(data=[go.Surface(z=z_data.values)])

fig.update_layout(title='Mt Bruno Elevation', autosize=False,
                  width=500, height=500,
                  margin=dict(l=65, r=50, b=65, t=90))

fig.show()

This code plots a 3-dimensional view on the example of CSV data. It can be scaled or rotated.

Any help?


texts: difference between two shapefile of one place

texts: I've two shapefile of one place. But they are not compatible. Actually when I add both of them into arcmap at the same time, there are difference between objects such as roads, parcels and etc. What is problem and how I can fix it?


texts: Running a non-QGIS Python 3 script within QGIS 3.4 processing toolbox just like that of QGIS 2.18

texts:
I've written around 30 scripts (both QGIS and non-QGIS scripts) in Python 2.7 and using it through QGIS 2.18 "Processing -> Scripts" tools so that within our office we can share and use the scripts.

We have moved to QGIS 3.4 and no issues with converting almost all of the 30 scripts from Python 2.7 to Python 3.7 and as of now, I am able to run those scripts through "QGIS Python Console". But i don't know how to convert those scripts to "Python processing tool scripts" which can be loaded in QGIS 3.4 and execute the scripts as I used to do in QGIS 2.18.

As an example, would like to brief about a script i am using - "Delete_lines_from_CSV"

1) Input will be a CSV file with 1st row - will have the list of CSVs for which the lines are to be deleted; IInd row - will have the corresponding new csv file name.

Below is the code I am using and it runs in the QGIS 3.4 Python Console without any issues.

import os, csv, re
from os import path, linesep
input_csvpath = r"D:\Test\pyqgis3\csvlist2.csv"
nlines = 2
nfirstlines = []

with open(input_csvpath, "r") as csvfile:
    reader = csv.reader(csvfile)

    for row in reader:
        input_csv= row[0]
        out_csv = row[1]

        with open(input_csv) as incsv, open(out_csv, "w") as outcsv:
            for x in range(int(nlines)):
                nfirstlines.append(next(incsv))
            for line in incsv:
                outcsv.write(line)

With the help of the links and comments mentioned in this thread, Writing a python processing script with Qgis 3.0 - I've modified the template code to the very basics

from PyQt5.QtCore import (QCoreApplication,
                      QVariant)
from qgis.core import (QgsProcessing,
                       QgsFields,
                       QgsField,
                       QgsProcessingException,
                       QgsProcessingAlgorithm,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterField,
                       QgsProcessingParameterFile,
                       QgsProcessingParameterString,
                       QgsProcessingParameterBoolean,
                       QgsProcessingParameterNumber,
                       )

import os, csv, re
from os import path, linesep

class deletelines(QgsProcessingAlgorithm):
    Input = 'input_iist'
    No_of_rows = 'no_of_rows'

    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def deletelines(self):
        return deletelines()

    def name(self):
        return 'Deletelines'.lower()

    def displayName(self):
        return self.tr('Delete_lines_from_CSV_TXT')

    def group(self):
        return self.tr('My_Scripts')

    def groupId(self):
        return 'myscripts1'

    def shortHelpString(self):
        return self.tr('Delete Lines from csv')

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterFile(name=self.Input,description="Input CSV List",extension="csv(*.csv)"))
        self.addParameter(QgsProcessingParameterNumber(self.No_of_rows,self.tr('Field length'), minValue=0, defaultValue=2))

and when I run this, I am getting

Don't know how to proceed further. In QGIS 2.18, I'd just paste the remaining code - as shown in the image below, highlighted in red and will get the outputs.

If I can learn how to implement this non-QGIS script within QGIS 3.4, can manage with the other non-QGIS, QGIS scripts.


texts: R Shiny Openstreetmaps - lock map within a set of upper and lower bounds

texts:
I am using Openstreetmaps with R/Shiny. I am trying to restrict map display to a given range of coordinates, i.e.xmin, ymin, xmax, ymax so that the users of the map won't be able to look other parts of the world other than the ones that are within the range.

I've made some research about it, but it hasn't been possible to find an answer. How can I do that?


texts: OpenLayers2 multiple maps -- which is used

texts:
My page contains multiple OpenLayers map objects loaded in the JS var domWidgetMaps . Each of them with a specific navigation.

var domWidgetMaps = $('div[id^="' + this.id + '"]');
for (var i = domWidgetMaps.length; i > 0; i--) {
    var options = {
        projection: mercator, //EPSG:3785/900913
        displayProjection: geographic, //EPSG:4326
        theme: null,
        /*
            Lines removed from OL 2.12
            maxExtent: new OpenLayers.Bounds(-20037508.34, -20037508.34, 20037508.34, 20037508.34),
            maxResolution: 156543.0339, //http://wiki.osgeo.org/wiki/WMS_Tiling_Client_Recommendation
        */
        units: 'm',
        controls: [
            new OpenLayers.Control.Navigation(),
            new OpenLayers.Control.PanZoom(),
            new OpenLayers.Control.ArgParser(),
            new OpenLayers.Control.Attribution()
        ],
        numZoomLevels: MAX_ZOOM_LEVEL,
        allOverlays: true
    };
    widgetMaps[i] = new OpenLayers.Map(this.id + i, options);
    var osm = new OpenLayers.Layer.OSM('osmBackground' + i, '', {visibility: true}, {buffer: 0});
    widgetMaps[i].addLayers([osm]);
    widgetMaps[i].zoomToMaxExtent();
}

If the user wants to compare several maps, features like: "zoom-in, zoom-out or pan" should be "shared".
The easiest solution (I think) is to identify which map I'm working on (How?) and to launch

var bboxSelected = widgetMaps[<SELECTED>].getExtent();
widgetMaps[2].zoomToExtent(bboxSelected);
widgetMaps[5].zoomToExtent(bboxSelected);
widgetMaps[6].zoomToExtent(bboxSelected);
widgetMaps[9].zoomToExtent(bboxSelected);

I tried

$('.map').click(function (e) {
});

It works if I click on the map, but it doesn't work if I click on zoom-in, zoom-out or pan buttons. OpenLayers intercepts the events.

How can I get the selected map based on zoom-in/out, pan?


texts: Open 4d NetCDF subdatasets with GDAL in Python

texts:
I have NetCDF files containing subdatasets with 4 dimensions: time, height, latitude, longitude.

Here is an example of the output of GDAL GetSubDatasets():

('NETCDF:"Martinique_ALT22.nc":height', '[30x30x30] height (32-bit floating-point)'),
('NETCDF:"Martinique_ALT22.nc":latitude', '[30x30] latitude (32-bit floating-point)'),
('NETCDF:"Martinique_ALT22.nc":longitude', '[30x30] longitude (32-bit floating-point)'),
('NETCDF:"Martinique_ALT22.nc":Water_Vapor_Concentration', '[2x30x30x30] Water_Vapor_Concentration (32-bit floating-point)')

When opening these subdatasets in Python with gdal.Open() and ReadAsArray(), the first two dimensions are overlapped and I get a 3 dimensions numpy array.

>>> band = gdal.Open(dataset.GetSubDatasets()[-1][0])
>>> array=band.ReadAsArray()
>>> print(array.shape)
(60, 30, 30)

I read somewhere that this is due to the GDAL raster format accepting only 3 dimensions: bands, rows, columns.

Is there a way to keep the first two dimensions separated and extract these subdatasets in 4-dimension numpy arrays?


texts:
This is more a theoretical question and not an actual problem that I have.

Is it possible to identify or classify projections by their numerical coordinates? Or alternatively rule out in some sort of a decision tree which projection some coordinates have? UTM and Lat-Long projection are normally easy to identify. What about the others?
Of all the possible (and likely many more) projection forms: http://en.wikipedia.org/wiki/List_of_map_projections

Assume the situation that you got a shapefile containing points without any .prj file or other information from the data provider. Is there some sort of pattern I could follow to rule out certain projections from the start?

Maybe someone can give some insight into the diversity and properties of coordinate projections.


texts:
I have a shapefile with a number of weed control blocks as polygons. They each have a unique block id number.
I am trying to join a non spatial csv file which has multiple entries for each block over 10yr period.
I have tried using Make a Query table but instead of coming up with 411 entries it comes up with over 19,000

I am using ArcGIS 10.1 for Desktop.


texts: Move files from Table Of Contents to a geodatabase in ArcPy?

texts:
I have a program that selects certain files and imports them to the Table Of Contents in ArcMap.

I would like to do some self intersect operations to all and put them both the original files and the intersected to a geodatabase.

Code attempt:

The below works the part after importing them to the Table Of Contents.
Can you help make the code functional?

for files in TOC: #somehow store all files located in Table Of Contents to a variable.
    arcpy.Intersect_analysis([[i, i],i + '_inte')#add the '_inte' to the output
    #and here move to the geodatabase both the original files and the intersected ones

My code to import to TOC:

arcpy.env.workspace = r"C:\Users\user\Desktop\20181018_CHECKS\05052"

walk = arcpy.da.Walk(datatype="FeatureClass")

for dir_path, dir_names, file_names in walk:
    for filename in file_names:
        if filename.endswith('.shp'):
            print(filename)

walk = arcpy.da.Walk(datatype="FeatureClass")
for dir_path, dir_names, file_names in walk:
    for filename in file_names:
        if filename.endswith('.shp'):
            arcpy.env.addOutputsToMap = True
            arcpy.env.overwriteOutput = True
            arcpy.MakeFeatureLayer_management(os.path.join(dir_path, filename),filename)

What is there now should be self-intersected and copied to the the geodatabase with the original files too.


texts: PostgreSQL - Creating DB, User privileges

texts:
I created a database named "wchiba" in postgres using  :

create database wchiba;

Then I create a user :

create a user saifeddine with password "0000";

Then I grant privileges :

grant all on database wchiba to saifeddine;

Everything now is ok

But when I load a shapefile to my table I have this  problem of password

FATAL:  password authentication failed for user "wchiba"

Also my database is wchiba not the user while the user created is saifeddine.

What's going on?

Does anyone know how to configure pg_hba.config


texts: Insert missing values in categorized style

texts:
In QGIS Bonn, when styling my layer to use the categorize option, I want to have all values within a range, for instance, 0 to 50000.

I am using the data from a CSV that has information about population density and I am using a grayscale to categorize it.

But the values are so different that 0 is white, 2000 is medium gray and 50000 is black.

This does not reflect the true differences between population per area and I do not want to add the missing integer values by hand, how do I do it?


texts: Does intersect result in non-simple geometry?

texts:
I have some features retrieved from SDE where I perform an Intersect on. Both feature-classes from which the features come from are polygon-featureclasses.

IFeatureCursor cur = myFeatureClass.Search(spFilter, false);
IFeature row = null;
while ((row = cur.NextFeature()) != null)
{
    ITopologicalOperator topOp = row.Shape as ITopologicalOperator;
    IArea baseGeom = (IArea)topOp.Intersect(feature.GetEsriGeometry(), row.Shape.Dimension);

    if (baseGeom is IGeometryCollection)
    {
        IPolygon4 mPoly = (IPolygon4)baseGeom;
        IGeometryCollection mGeomCol = null;
        try
        {
            // this may throw a COM-exception when the polygon is not simple (errorcode 0x80040218 = E_GEOMETRY_NOTSIMPLE)
            mGeomCol = (IGeometryCollection)mPoly.ConnectedComponentBag;
        }
        catch
        {
            if (baseGeom is IGeometryCollection) mGeomCol = (IGeometryCollection)baseGeom;
        }
    }

    // do something with our collection
}

As of the ESRI-docs we know that performing an Operation on ITopologicalOperator will always return a topological correct geometry. Therefor both the IsKnownSimple- and the IsSimple-properties for baseGeom are true.

When I call ConnectedComponentBag from IPolygon4 it throws a COM-exception:

  (0x80040218) The operation cannot be performed on a non-simple geometry.

That makes me wonder how the property of a topological correct geometry may become topological incorrect. The docs for the property state that it consists of the rings of the original polygon, so where is the difference between my baseGeom which is topological correct and my ConnectedComponentBag which is not?


texts: Find the convex hull around a set of polygons?

texts:
I'm using Postgres 9.4 with PostGIS 2.1. I have a table with the following columns:

    Column    |          Type
--------------+------------------------
 ogc_fid      | integer
 wkb_geometry | geometry(Polygon,4326)
 plot_id      | integer
 org          | character varying

I want to get all the polygons relating to a particular organisation, draw the convex hull around them, then calculate the area of that convex hull.

I've tried this query:

SELECT
 SUM(ST_AREA(ST_CONCAVEHULL(wkb_geometry, 0.1)::geography)) AS area
 FROM mytable
 WHERE org='myorg'

But it seems to be first finding the convex hull of each individual polygon, then adding together the areas. This is not what I want.

How can I edit this query to first find the convex hull of the entire polygon set, and only then calculate the area?


texts: QGIS or Postgis - How to split polylines with a point layer?

texts:
In QGIS, I would like to split the polylines in a layer loaded in the workspace with the points in another loaded layer. If a point intersects a polyline, this polyline should be transformed into two polylines...

I've read that this could be done with grass but i have no idea how to do it with layers (SHP)loaded in the QGIS workspace.


texts:
I have a large project with provincial data. The project is currently using a data driven map with a clipped data frame for the province. The issue is the data I'm trying to show is thousands of point labels within each province, and the now clipped border is clipping labels as well. How do I stop the labels from being clipped along with the data frame?

Also, I only want to show labels that are applicable to that specific province and not the neighbouring provinces.


texts: Why is arcpy.CalculateField_management() writing 1899-12-30 00:00:00 instead of expected date?

texts:
I'm using ArcGIS 10.3.1 for Desktop to test before writing an answer to Entering date with arcpy.GetParameterAsText()? I have written the test code below that I am running against a test feature class (testFC) with four line features to try and update its dateField with a constant date value of "01.07.2015".

import arcpy

dmyString = "01.07.2015"
fc = r"C:\temp\test.gdb\testFC"
item = "dateField"

dmyDate = datetime.datetime.strptime(dmyString,'%d.%m.%Y')
print dmyDate

arcpy.CalculateField_management(fc,item,dmyDate,"PYTHON")

with arcpy.da.SearchCursor(fc,[item]) as cursor:
    for row in cursor:
        print row[0]

I think the above should work, and from the output below you can see that the date object (dmyDate) seems to be created fine, but the values written to each of the four features are not what I would expect.

>>> ================================ RESTART ================================
>>>
2015-07-01 00:00:00
1899-12-30 00:00:00
1899-12-30 00:00:00
1899-12-30 00:00:00
1899-12-30 00:00:00
>>>

Can someone explain why 1899-12-30 00:00:00 is being written instead of 2015-07-01 00:00:00?

The attribute table displays as below:

If I replace the CalculateField line above with either:

arcpy.CalculateField_management(fc,item,"date '2015-07-01 00:00:00'","PYTHON")

or

arcpy.CalculateField_management(fc,item,"date '2015/07/01 00:00:00'","PYTHON")

as per @Vince's comments I receive an error:

Traceback (most recent call last):
  File "C:\temp\test.py", line 11, in <module>
    arcpy.CalculateField_management(fc,item,"date '2015-07-01 00:00:00'","PYTHON")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\management.py", line 3457, in CalculateField
    raise e
ExecuteError: ERROR 000539: SyntaxError: invalid syntax (<expression>, line 1)
Failed to execute (CalculateField).


texts:
How to make a halo around the line of one pixel using OpenLayers?

For example, the line (vector) has  a blue color and its frame (one pixel) has a black color.


texts: OpenLayers deploy: what libs does selectFeature rely on?

texts:
I am trying to minimize my code by modifying OpenLayers deploy configuration file in order to reduce its load time. However, It seems that I can't include all necessary libs that OpenLayers.SelectFeature.Control needs. Here is the config file:

    [include]
OpenLayers/Map.js
OpenLayers/Kinetic.js
OpenLayers/Projection.js
OpenLayers/Handler.js
OpenLayers/Control.js
OpenLayers/BaseTypes/Pixel.js
OpenLayers/BaseTypes/LonLat.js
OpenLayers/BaseTypes/Bounds.js
OpenLayers/Layer/SphericalMercator.js
OpenLayers/Layer/XYZ.js
OpenLayers/Layer/WMS.js
OpenLayers/Layer/Vector.js
OpenLayers/Layer/ArcGISCache.js
OpenLayers/Control/TouchNavigation.js
OpenLayers/Control/Geolocate.js
OpenLayers/Control/Attribution.js
OpenLayers/Control/SelectFeature.js
OpenLayers/Control/DrawFeature.js
OpenLayers/Control/ModifyFeature.js
OpenLayers/Control/Snapping.js
OpenLayers/Control/Split.js
OpenLayers/Control/Panel.js
OpenLayers/Control/Button.js
OpenLayers/Control/PinchZoom.js
OpenLayers/Feature/Vector.js
OpenLayers/Feature/WFS.js
OpenLayers/Handler/Click.js
OpenLayers/Handler/Feature.js
OpenLayers/Handler/Hover.js
OpenLayers/Handler/Box.js
OpenLayers/Renderer/SVG.js
OpenLayers/Renderer/Canvas.js
OpenLayers/Format/GeoJSON.js
OpenLayers/Format/KML.js
OpenLayers/Format/WFS.js
OpenLayers/Format/WFSCapabilities.js
OpenLayers/Format/WFSCapabilities/v1.js
OpenLayers/Format/WFSCapabilities/v1_0_0.js
OpenLayers/Format/WFSCapabilities/v1_1_0.js
OpenLayers/Format/WFSDescribeFeatureType.js
OpenLayers/Format/XML.js
OpenLayers/Request/XMLHttpRequest.js
OpenLayers/Protocol/HTTP.js
OpenLayers/Protocol/WFS.js
OpenLayers/Protocol/WFS/v1.js
OpenLayers/Protocol/WFS/v1_0_0.js
OpenLayers/Protocol/WFS/v1_1_0.js
OpenLayers/Strategy/Fixed.js
OpenLayers/Strategy/BBOX.js

When I test selectFeature feature which has a feature info popup in OnSelect, it always tells me following error:

cannot call "add" of undefined

which I think it points out that in the file of selectFeature.js, under "pixelToBounds", the variable "pixel" is undefined, and then I assume it's relevant to the recognition of "evt".

could somebody give me some hints?

thank you!


texts: Need to understand desc.workspaceType to use in conditionals

texts:
I need to differentiate between two possible types of workspaces for an arcpy script. So I check the workspaceType:

rootWorkspace = raw_input('Geodatabase workspace to analyze: ')
desc = arcpy.Describe(rootWorkspace)
wsType = desc.workspaceType
print 'Workspace Type: {}'.format(wsType)

which seems to work fine and prints the correct type. But, if I try to then check the type (so I can, eventually, do different things depending on what it is):

if wsType is 'FileSystem': print 'file system'
if wsType is 'LocalDatabase': print 'gdb'
if wsType is 'RemoteDatabase': print 'sde'
else: print 'not reading this correctly'

I always get not reading this correctly (whether I'm inputting a directory or geodatabase).

According to Arc help about Workspace properties, workspaceType should be a string result (not an Object etc.), so I am confused why a string comparison isn't working.


texts: Will killing ArcMap process while it's copying features corrupt the file geodatabase?

texts: So..., I made the unfortunate decision to try and copy/paste 22,000+ address points into the same file geodatabase.  It's taking soooooooo long....  So my question is, how much of a chance is there to corrupt a file geodatabase during a copy/paste operation if you kill the ArcMap process midway thru?


texts: NetCDF hourly data to year raster layers (using Python)

texts:
Having a NetCDF with a 'time', 'lat' and 'lon' dimension, where time is hourly during a year and the variable to extract is temperature, how do I calculate mean temperature throughout the year using Python?
(using Python and netcdf4 or xarray or numpy)


texts: OpenLayers ECQL BBOX and INTERSECTS: Can I set SRS name?

texts:
I am using GeoServer and am trying to understand one problem:
Using this statement returns successfully the polygon intersecting with my BBOX. The layer is EPSG:2193

http://myserver/fwsys/wms?service=WFS&version=1.0&request=GetFeature&typeName=fwsys:fwsys_region&CQL_FILTER=BBOX(the_geom,19438298.842,%20-4502716.974,19702465.211,%20-4587714.950,%27EPSG:3857%27)

What I would like to point out here is that the layer is in EPSG:2931 - my request is in EPSG:3875. This works well

Now I followed the instructions and tried to do the same for a single point:

http://myserver/fwsys/wms?service=WFS&version=1.0&request=GetFeature&typeName=fwsys:fwsys_region&CQL_FILTER=POINT(the_geom,19438298.842,%20-4502716.974,%27EPSG:3857%27)

This returns the following error output:

<ServiceExceptionReport version="1.2.0"
  xsi:schemaLocation="http://www.opengis.net/ogc http://schemas.opengis.net/wfs/1.0.0/OGC-exception.xsd">
  <ServiceException>
    Could not parse CQL filter list.
    Encountered "POINT ( the_geom" at line 1, column 1.
    Was expecting one of:
    "not" ...
    "include" ...
    "exclude" ...
    "(" ...
    "[" ...
    "id" ...
    "in" ...
    <IDENTIFIER> ...
    <DATE_TIME> ...
    <DATE> ...
    "-" ...
    <INTEGER_LITERAL> ...
    <FLOATING_LITERAL> ...
    <STRING_LITERAL> ...
    "true" ...
    "false" ...
    "point" ...
    "point" "(" "-" ...
    "point" "(" <INTEGER_LITERAL> ...
    "point" "(" <FLOATING_LITERAL> ...
    "point" "not" ...
    "point" "like" ...
    "point" "ilike" ...
    "point" "exists" ...
    "point" "does-not-exist" ...
    "point" "is" ...
    "point" "between" ...
    "point" "=" ...
    "point" ">" ...
    "point" "<" ...
    "point" ">=" ...
    "point" "<=" ...
    "point" "<>" ...
     Parsing : POINT(the_geom,19438298.842, -4502716.974,'EPSG:3857').
  </ServiceException>
</ServiceExceptionReport>

Is it possible to use 'POINT' in the way I do? What am I doing wrong?


texts: How to clip a 3-band raster when unwanted, bordering cells have values?

texts:
I have myself a large raster (ecw) that comprises of 3 bands. I'm trying to find the most convenient way of clipping away unwanted black areas around the actual imagery itself. See the raster in question below.

The problem is that surrounding black cells are not classified as nodata, but have differing values such as 0/0/0, 1/0/1, 2/0/1...Thus it takes take time and effort to choose all of them in e.g. raster calculator. Also these values might exist in the area of interest.

For example, if I would set value 1/0/1 to be transparent in layer properties --> transparency dialog, I get the result shown below.

Does anyone have any tips about how to proceed with this? Am I forced to manually create a masking polygon?


texts: Which the columns available on the TIGER geocoder

texts:
I have a local setup of the tiger geocoder and I want to extract as much information as possible from it.
Currently, from a query

SELECT g.rating, ST_X(geomout) As lon, ST_Y(geomout) As lat, (addy).* FROM geocode('75 State Street, Boston MA 02109', 1) As g;

I'm only getting this features:

rating |        lon        |       lat        | address | predirabbrev | streetname | streettypeabbrev | postdirabbrev | internal | location | stateabbrev |  zip  | parsed
--------+-------------------+------------------+---------+--------------+------------+------------------+---------------+----------+----------+-------------+-------+--------
      0 | -71.0556722990239 | 42.3589914927049 |      75 |              | State      | St               |               |          | Boston   | MA          | 02109 | t

Which other information is available to extract?
I'm looking for census track, building type, parcel number...


texts: A basic question on layers, sublayers, and cartodb.js

texts:
I am really new to cartodb and have spent a number of hours trying to figure out the distinction between layers and sublayers.  Specifically within cartodb.js.   For example, I have a map with two layers that I created using cartodb.  I then use cartodb.createVis to get at it from javascript.   Most things work, but to set an infowindow on a sublayer, I seem to have to do:

layers[1].layers[1].sub.infowindow.set(stuff)

but, when I define the infowindow's template, only some of the template tags seem to be populated.  Any ideas?  In particular:

1)  Is there clear documentation somewhere on layers, sublayers, events, and infowindows somewhere?

2)  Is there any way within an infowindow to inspect {{content.data}} to find out what fields it actually has?


texts: Static maps API that serves roads/labels without basemap?

texts:
I am looking for a static maps API (similar to https://www.mapbox.com/developers/api/static/ ) that can serve the map without an underlying base layer. My intention is to generate images with the roads/labels layer on top of a custom base layer.
Does such a thing exist? If not, what can I use to create it?


texts: Using satellite background images to draw features in PostGIS database?

texts:
What tool do I need to have satellite background images to draw features of a PostGIS database?

I looked about Geomajas but seems very complex for what I'm trying to do.

I have osm data imported in a postgis database, as a layer in QGIS. The second layer is what I am drawing, and the third layer that I ask here is how to have satellite imaginery on bottom of the two mentioned ones as a visual aid in drawing features.


texts: Configurable Map Viewer (CMV) - Open Help on page load

texts:
I'm wondering, as opposed to opening the help dialog when the help link is clicked:

var help = domConstruct.place(this.html, this.domTarget);
on(help, 'click', lang.hitch(this.parentWidget, 'show'));

How might I go about loading the help dialog when the viewer first loads?  So, by default the help opens and needs to be actively closed to interact with the viewer?

It looks like I could use this.parentWidget.show(); in postCreate but it doesn't block the ability to interact with the underlying viewer.

How is that accomplished?


texts: Transfer RPC metadata from one geotiff to another

texts: I am performing an image processing task on one geotiff that produces another, however the metadata, including the RPCs, aren't propagated from the input to the output. Is there a way to copy the RPCs from the one image to the other, using GDAL or similar?


texts: If vb code is also possible I could give that a try aswell.


texts: Selecting in expression features with line break

texts:
I'm trying to filter features in a QGIS layer to those with the following in their "description" attribute.
CC_W

3

I think what I need is something like this, with the appropriate escape sequences replacing LINEBREAK,LINEBREAK:
"description" LIKE '%CC_W	LINEBREAK,LINEBREAK3%'
I've tried \r and \n to no avail. Both of them generate no filter hits even when I search them alone over my feature set, which definitely does have line breaks in the "description" attributes.
"description" LIKE '%\r%' 0 hits
"description" LIKE '%\n%' 0 hits


texts: Distance vs cost in QGIS network analysis

texts:
I'd like to calculate shortest route and distance in QGIS 3.2. In the past I have used the v.net.iso, v.net.distance and v.net.path tools. However, I've noticed that there is no distance parameter there, but cost. If I want to calculate the distance in meters where do I specify this?
Do cost mean distance on the screenshot? How do I set it to distance?


texts: Missing results from ReshapeGeometry method on PyQGIS

texts:
I am trying to use the method reshapeGeometry from QgsGeometry class in a script to reshape a feature geometry (polygon) according to a line.
The result of the process is 0, so positive, but I can't see any change in my geometry.
This is an example of my script:
layer = iface.activeLayer()
line = QgsLineString([QgsPoint(1, 1), QgsPoint(2, 2)])

for feat in layer.selectedFeatures():
    c = feat.geometry().reshapeGeometry(line)
    print(c)

The feat in this script is a simple polygon with WKT.
Polygon ((1.16531548008501851 1.44378001491727948, 1.60296509683766097 1.93729979551068476, 2.06544589751953112 1.57103983255457269, 1.60296509683766097 1.0309615820938649, 1.16531548008501851 1.44378001491727948))
0

as observable in the image below (in green). The red line is the representation of the line used for the reshaping.

I have tried both to make the layer editable before and saving the edits after and to reload the layer after the operation.
I am working in the PyQGIS console on QGIS 3.16 and 3.19 dev.


texts:
I am attempting to load static tiles from this ArcGIS MapServer into R, and generate static maps with tmap 3.3.
The package has a pretty convenient way to do this in view-mode, generating interactive maps.
url_tile_services <- "https://tiles.arcgis.com/tiles/RTK5Unh1Z71JKIiR/arcgis/rest/services/"
tile_example <- "Global_Dam_Watch/MapServer/tile/{z}/{y}/{x}"

data("World")

## Interactive map
tmap_mode("view")

tm_shape(World) +
  tm_tiles(paste0(url_tile_services, tile_example), group = "LABELS") +
  tm_symbols(size = "gdp_cap_est")

Conversely, the method to map static tiles in plot-mode, i.e. returning static maps, is rather unclear. In fact the code below returns the symbols layer but not the tiles layer. Unfortunately, neither a Warning nor an Error message is returned, and this makes the troubleshooting more challenging because I have no idea of what's going wrong.
url_tile_services <- "https://tiles.arcgis.com/tiles/RTK5Unh1Z71JKIiR/arcgis/rest/services/"
tile_example <- "Global_Dam_Watch/MapServer/tile/{z}/{y}/{x}"

data("World")

## Static map
tmap_mode("plot")

tm_shape(World) +
  tm_tiles(paste0(url_tile_services, tile_example), group = "LABELS") +
  tm_symbols(size = "gdp_cap_est")

This question and this issue raised on the GitHub deal with this matter but all suggested solutions focus on OpenStreetMap tiles, while I need several from those provided in the MapServer linked above to produce static maps.
Any idea how to accomplish that?


texts: Define a default style to workspace in GeoServer

texts: I would like to know if there's a way to define a default style to a whole workspace in GeoServer? For instance each time a raster is published will have the default style we defined using YSLD.


texts: How to remove a QGIS project from the "Recent Projects" list?

texts:
I have a project in the QGIS Recent Projects list that doesn't exist anymore, but it keeps being listed (greyed out and without thumbnail image). Is there any way to remove it from the list? Perhaps by manually editing a history file or something like that?

I've found a feature request for this, and a related discussion, but I'm wondering if there's anything I can do now in QGIS 2.14.


texts: Trouble with GeoJSON in HTTP GET Request

texts:
I using a Vector Layer with the HTTP protocol to request GeoJSON from my web server, and am trying to get this GeoJSON to display on my map. The GET request appears to be working correctly, as I can see the example JSON I requested in Firebug. However, I am getting an error telling me that "this.format is null     Openlayers.js (line 1186)". Unfortunately, the OpenLayers.js file is essentially unreadable so I am having a hard time figuring out what my problem might be.

I set the format for the Vector Layer to GeoJSON. My code for the layer is below:

var newLayer = new OpenLayers.Layer.Vector("Filtered by Zoom Level", {
                        strategies: [new OpenLayers.Strategy.BBOX(),
                                    refresh,
                                    ],

                        protocol: new OpenLayers.Protocol.HTTP({
                            url: "http://server:8080/test-webapp/servlet",
                        }),

                        filter: filter,

                        styleMap: styleMap,
                        format: new OpenLayers.Format.GeoJSON({
                            externalProjection: new OpenLayers.Projection("EPSG:4326"),
                            internalProject: new OpenLayers.Projection("EPSG:900913")
                        })

                    });

Does anybody have any ideas as to what the issue might be?

Thanks,

-Rob


texts: How to use consistent symbology between multiple stretched rasters in ArcMap?

texts:
I have a series of continuous rasters representing the same phenomenon over time, so I need to use the same symbology across them all.

So I set the stretched symbology of the first raster, with a color ramp from the min/max values:

Now in the second raster, I import the symbology from the first raster. Unfortunately, although the labels match, the actual min/max values are not the same as the first layer:

So the red in the second raster corresponds to a value of 0.47, rather than 0.66.

How can I ensure that both rasters are using a consistent symbology?


texts: Divide broken line layer in seperate clusters based on location

texts:
Building on a earlier asked question on this forum I would like to ask an additional question in getting the task done I asked here: Dividing polyline layer in clusters with specific requirements in ArcGIS for Desktop?.

To do some preparation for this task I would like create a script to divide my line layer in separated parts making a cluster if the links are connected. So, in the next picture I show an example of what I mean.

To get this done I created a spatial join 'one to many' and match option 'boundery touches'. So that creates a table with all values that are connected to each other. My Idea is to take the FID of the 'original table' (for example 0), and search that value in the 'join table' in the field 'Target_FID'. Then take the values from the JOIN_FID column in the same rows (47, 109, 111). And then search these values (47, 109, 111) again in the Target_FID' column and get the values from the JOIN_FID column in the same rows (not visible in the picture). And do this until I get a list of all the values (FIDs) which are connected to each other. Next step then is to assign the same value in a new column to all the FIDs in this list in the 'original table'. To visualize this:

If the list of FIDs which are connected to each other was for example (0,1,2,3,4,5 and 6) give this value 1. And so on for other list of FIDs which are connected giving them values 2,3,4 etc....
This should create clusters of the connected parts.

To get this done was thinking of creating a script using searchcursor and updatecursor. I searched on how to get it done with these options but I couldn't figure out how to do it for my typical situation. Also because I am new to using these options in scripts. So, I am wondering if using the search and update cursors is the right method to get this task done? And if yes does anyone have a idea how the python script should look like? Or if it is not the right method, does anyone have a better idea to get this task done? It should not be to difficult I would say, but just haven't found the right method yet.


texts:
I am familiar with XY, XYZ, and XYZM, features but I created a spatial layer in PostgreSQL and added it to my QGIS. When I opened the Vertex Editor it showed me values for XYR and I am not familiar with the R; can anyone tell me what it is? Is it the equivalent of ESRI's M value?
I looked for a QGIS forum but their website directed me here.


texts: What are the units of velocity in Google Location History JSONs?

texts:
I have downloaded and used different code to visualize the data. Some of the code works with speed as it is available in the JSON files downloaded using Google Takeout. The speeds seem to be too slow to be in km per hour.

I have found no documentation on the topic. I filed a google help ticket with no answer. Aside from this, everything works great, but I need this information to be able to interpret my results.


texts: Merging lines with slightly different endpoints in QGIS

texts:
In QGIS, I have two lines. The endpoints are near each other, but not quite aligned. I'd like to turn the two lines into one, with the two endpoints merged into midpoints.

The closest I've been able to come is with "Merge selected features." That gives me one feature, but it's comprised of two disjoint lines -- there's no line between the endpoints.

I'm looking for something that behaves similarly to Illustrator's Join Paths feature (Object > Paths > Join).

Does QGIS do that?


texts: Seeking equivalent of ArcGIS Pro Summarize Within in QGIS

texts:
I have two polygon layers : "Building foortprints" and "Building buffer" (which is
a buffer around each building and they do overlap each other). I want to calculate the sum area of all buildings in each buffer. I tried to use Summarized Within tool in ArcGIS Pro, but it gives me error 99999. I guess because the buffers overlaps, it cannot run the tool succesfully.

Is there any way that I can do this calculation in QGIS?

For example in the below image, for each building's buffer I should have the sum of three other buildings area.


texts: PostGIS Voronoi Polygons with extend_to parameter

texts:
I am creating voronoi polygons for locations inside a country, using PostGIS.
The function creates these areas for ~80,000 locations in less than a second.
However, I can only get the polygons to extend to the default envelope area, which is 50% of the bounding box of points.
WITH voronoi (vor) AS
 (SELECT ST_Dump(ST_VoronoiPolygons(ST_Collect(geom))) FROM public.branches_atms_0601)
SELECT (vor).path, (vor).geom FROM voronoi;

This generates the polygon(s) below:

What I would like it to do is to clip this to the footprint of UK. This can apparently be defined within the function https://postgis.net/docs/ST_VoronoiPolygons.html
However when I try to run it with my UK footprint polygon, I get the error with the code...
WITH voronoi (vor) AS
 (SELECT ST_Dump(ST_VoronoiPolygons(ST_Collect(geom), public."Dissolved")) FROM
public.branches_atms_0601)
SELECT (vor).path, (vor).geom FROM voronoi;

ERROR:  missing FROM-clause entry for table "public"
LINE 2: ...LECT ST_Dump(ST_VoronoiPolygons(ST_Collect(geom), public."Di...
                                                         ^
SQL state: 42P01
Character: 82

How can I fix this? I assume the issue is coming from assigning the extends_to parameter incorrectly, but I don't know how to approach it...


texts: I have installed the SRTM downloader, it calculates the number of tiles and then start the download in a flash, but without actually download. It does not ask to log in... and proceeds!


texts: Taking only relevant rows from database, based on common codes in QGIS?

texts:
I have a shp file that is missing some data, e.g. attributes. I have the missing data in a csv file. How can I get ONLY the relevant rows from the csv file added to the current shp file's data table? I.e. ONLY those who match the rows that are already in the shp file.

The good news is yes, they both have the same ID code in a column of each, so 'lining them up' should be easily do-able.

I am using QGIS.


texts: Having trouble starting an edit session on SDE with edit permissions set for user

texts:
I have at least 15 people editing on one of our Enterprise SDE Geodatabases. I set up another user to edit and it is giving me an error I can't figure out. I set the same exact permissions as the other users who are editing fine and have been for months.

Database is SDE on SQL Server 2008, Arcgis Server 10.3, he is using ArcMap 10.3.1

Sidenote he had just got Windows 10, I think everyone else is using Windows 7.


texts: Calculate potential of solar radiation on rooftops (without LiDAR data)

texts:
I'm trying to calculate the solar radiation for rooftops, but I have a big problem, I don't have LIDAR data. I need to give the buildings the information of height. I've tried to do these, but for a building I'm going to get for example three types of heights when I should only get one.

So, I want to ask if anyone could suggest another way to calculate the buildings height and the surface of rooftops, so that I can calculate the potential of solar radiation for rooftops without using LiDAR data.


texts: How to label using just numbers from start of text field?

texts:
In ArcMAP 10.2, I am trying to format the labels on my map.

If the label is

1402 Bowie Road

But  I want only to display numbers:

1402

What is the expression to create this?


texts: What is FME equivalent to ArcGIS Dissolve?

texts: I am trying to dissolve polygons that are spread out across a region. ArcMap's ModelBuilder will dissolve these polygons into a single, multi-part feature. How can I get FME to do this same thing? I have reviewed the Dissolver and FeatureMerger transformers without any luck.


texts: Create a point shapefile from table in .accdb database

texts:
I am a beginner in Python.

I have a database in .accdb format of Access 2007. The database has a table called "T_point_RT" which has point data.

I want to create a point shapefile from North and Est coordinates (in T_point_RT table) in a way that each record in the shapefile corresponds each point. The shapefile must have also the ID field for each point.

This is the script, the problem is that the shapefile that is created has only one record/one point (the last one). What is the problem ??

import arcpy
import pyodbc
import numpy
import datetime

data=datetime.date.today().strftime('%Y%m%d')

conn = pyodbc.connect(r'Driver={Microsoft Access Driver (*.mdb, *.accdb)};DBQ=C:\Python_Elaborazioni\Dati\Database\DB_COP_Field_RedTeam.accdb;')

cursor = conn.cursor()

for row in cursor.execute('''select ID_point, Est, North from T_point_RT'''):
    array = numpy.array([(row[0], (row[1], row[2]))],
                   numpy.dtype([('Id_Point',numpy.unicode_,6),('XY', 'f8',2)]))
    SR = arcpy.SpatialReference(32632)
    SHP_WGS84= "C:/Python_Elaborazioni/Dati/Shp/P_DB_COP_Field_RedTeam_"+ data +" _wgs.shp"
arcpy.da.NumPyArrayToFeatureClass(array, SHP_WGS84, ['XY'], SR)

outFC_MM = "C:/Python_Elaborazioni/Dati/Shp/P_DB_COP_Field_RedTeam_"+ data +"_gb.shp"
SR_MM = arcpy.SpatialReference('Monte Mario Italy 1')
arcpy.Project_management(SHP_WGS84,outFC_MM,SR_MM)
print "Finito"


texts: Can image layer be rotated in OpenLayers?

texts:
While working on Ground Overlays, I noticed that it has a tag called <rotation>

 <?xml version="1.0" encoding="UTF-8"?>
    <kml xmlns="http://www.opengis.net/kml/2.2">
      <GroundOverlay>
        <name>Large-scale overlay on terrain</name>
        <visibility>1</visibility>
        <description>Overlay shows Mount Etna erupting on July 13th, 2001.</description>
        <Icon>
          <href>http://kml-samples.googlecode.com/svn/trunk/resources/etna.jpg</href>
        </Icon>
        <LatLonBox>
          <north>37.91904192681665</north>
          <south>37.46543388598137</south>
          <east>15.35832653742206</east>
          <west>14.60128369746704</west>
         <rotation>11.9556640799496235</rotation>
        </LatLonBox>
      </GroundOverlay>
    </kml>

My question is - can an image layer be rotated in 2D view?  Following is the code I am using to add an Image and set its bounding box.

var options = {
            'opacity': 1.0,
            'isBaseLayer': false,
            numZoomLevels : 20

        };
        var image = new OpenLayers.Layer.Image(
            groundOverlay.name,
            groundOverlay.href,
            new OpenLayers.Bounds( groundOverlay.west, groundOverlay.south,  groundOverlay.east, groundOverlay.north),
            new OpenLayers.Size(1, 1),
            options
            );

        map.addLayers([image]);

Following the screen shot when above KML is displayed on my local system without rotation.

Google Maps does not seem to rotate the image, but Google Earth plugin does.


texts: mile distance around each feature ArcGIS

texts:
I have a school layer information and a route map layer. In ArcMap, how can I show the x mile distance around each school base on route information.

I did some research and I found different approaches such as Network Analysis tool, Proximity Tool, Buffer Geo-processing tool,... .

I am not sure what approach I need to use.


texts:
It is said that unit of measure for EPSG:3857 is the metre.

How it can be, if the range of these coordinates is

Projected bounds:
-20,026,376.39 -20,048,966.10
20,026,376.39 20,048,966.10

(from the same source).

This means approximately a square of 40,000 x 40,000 kilometers. This length is correct for the equator, but how can it be in general if the circumference at high latitudes is smaller?

UPDATE

The length of the parallel at latitude 85.06 is

$R_earth * cos(85.06) * 2 * \pi$ = 3454 kilometers and it is << 40000

My question is: is this parallel projected onto the entire upper edge of the 40000 x 40000 square or onto part of it 3453m long?

In the first case the unit of measure IS NOT METRES, since it varies from metres at the origin to 3453/40000 = 9 centimeters at parallel 85.05.


texts: PermissionError : [WinError 32] --> how to close raster file handler after use?

texts:
I created a function where I access to a Raster folder and get all files CRS

            for Rsubdir, Rdirs, Rfiles in os.walk(Rrootdir):
                for Rfile in Rfiles:
                    Rfile, Rextension = os.path.splitext(os.path.join(Rsubdir, Rfile).replace("\\", "/"))
                    if Rextension == '.tif':
                        RfileName = os.path.join(Rsubdir, Rfile).replace("\\", "/") + ".tif"
                        if os.path.isfile(RfileName):
                            RfileInfo = QFileInfo(RfileName)
                            RbaseName = RfileInfo.baseName()
                            Rlyr = QgsRasterLayer(RfileName,RbaseName)
                            coordR=Rlyr.crs().authid()
                            RasterCRS.append(coordR)

the problem is when I try to access the same folder files after this function call, I get a PermissionError

PermissionError: [WinError 32] Le processus ne peut pas accéder au fichier car ce fichier est utilisé par un autre processus: 'chemin_fichier.tif'

which means that the files are always opened. is there any way to close them  in the end of the function?


texts:
I am relatively new to QGIS.

I am looking at the percentage of an population that has geographical access to healthcare.
Therefore I have one shapefile with village boundaries including population number of each village.
Another shapefile is all Healthcare facilities in the area buffered with 5km zones.
I now want to calculate the percentage of the village area that is covered by the Healthcare facility buffers and multiply this by the my population numbers.

I did manage to create an intersect shapefile of the two formal shapefiles that created correct new polygons of intersection per village. I now, however, do not know how to sum the areas (Area_1) of all these buffer polygons which are coded equally in the attributable table (under P_06_ID) and divide this by the total village area (AREA). And then multiply this number by the population (Population).

The command does indeed calculate the percentage of the population that overlaps with that specific polygon. However it still does not sum the different polygons.

So what I now want to do is to sum those numbers with the same parish ID automatically (I do have around 500 parishes). So in this example sum PopPerc from cell 699- 709.

Any suggestions?


texts:
I have one single shapefile with many features. There are three different topological errors that I want to be identified (only identified).

Line A consists of many segments and some of segments overlap each other;
Line 2,3 and 5 don't touch/snap exactly line A as it should be;
Line 4 intersects line A and I want it split into two segments where it intersects with main line A.

Solution I tried:

"LineonLineOverlayer" and "Snapper" are working pretty good the only issue is that they both return not only overlaping respectively not touching segments but everything which means also the segments that do not have problems like 1 and 6! I have a bigger problem with "Intersector" which should be able to give me just line 4 which is not the case.

Any advice ?


texts: Feature Layer (Polygon) Rendering issue in ArcGIS Server 10.1

texts: I have a feature layer consisting of 96,000 features, but we are displaying a maximum of 13,000 features in our map application. There is an issue that for some zoom levels the features are visible and for other, intermediate zoom levels the features are not visible. I have deleted the feature class and restarted the service in case there is issue with the feature class, but that also did not help.Thank you in advance and any kind of help will be appreciated.


texts: Coordinate reference system transformation

texts:
I am looking to transform the CRS of the geometry from my initial shapefile data to standard lat-long coordinates - espg 4236.

Here's the geopandas dataframe:

I wrote some script below that transforms the geometry to the espg CRS.

geodata.crs

{'proj': 'lcc',
'lat_1': 37.06666666666667,
'lat_2': 38.43333333333333,
'lat_0': 36.5,
'lon_0': -120.5,
'x_0': 2000000,
'y_0': 500000.0000000002,
'ellps': 'GRS80',
'units': 'us-ft',
'no_defs': True}

pointsEB is point datatype:

POINT (6091769.359598021 2090180.343520682)

# Function to transform CRS projection
inProj = Proj(geodata.crs)
outProj = Proj(init='epsg:4326')

pointsEBt = []

for pt in pointsEB:

    x,y = transform(inProj, outProj, pt.x, pt.y)
    pointsEBt.append(Point(x,y))

Transformed Point:
pointsEBt[1].x, pointsEBt[1].y
(-68.42841348959915, 40.16978821064519) --> Located in Antartica.

My question: There seems to be something going on with the conversion to lat-long
 CRS. I used this as reference and "+init=epsg:4326" seems to be the format.
 However, when I plot a transformed point on the map the location is not what I
 expect it to be, far from it.


texts:
In QGIS 2.10.1 I want to create a point layer out of my polyline shapefile and all the points shall have a fixed distance (in meters) while I use WGS84 as coordinate system.
So what do I use for this? I tried QChainage but it does not work since my coordinate system is WGS84 and thus I have only the option to set it in degree.


texts: How to use PostGIS to handle complex geoprocessing workflows?

texts:
Our organization is considering moving our geoprocessing workflow to PostGIS. We are currently using ArcGIS, with a plethora of custom Python tools used in ModelBuilder. We're moving most of our data into PostGIS to be consumed by a variety of apps, and now we're asking if it also makes sense to perform the data processing there as well.

We process data to be compatible with our software. A customer purchases our software, gives us their data, and we process it to be optimized for use in our software. This requires us to build a variety of tools to handle varying qualities of input data. We can't expect to receive data in a particular format or schema, so we build tools to map input fields to output fields, parse single fields into multiple fields, merge multiple datasets, etc. We also perform spatial joins, intersections, trim whitespace and concatenate fields, and many other common operations. PostGIS appears to be perfectly capable of performing all of our processing needs.

For those of you who use PostGIS to do your data processing, do you have any advice for organization, tools to use, etc.?

do you use it in conjunction with QGIS python processing?
are people using a Python ORM for non-web processing? I've been leaning towards using GeoDjango since it has a Python ORM for PostGIS. Our initial test of using PostGIS to process data has many large SQL text blocks in Python code and we are thinking that the GeoDjango ORM may help with creating more manageable and readable code. There's also the GeoAlchemy ORM that interacts similarly with PostGIS, and doesn't appear to be as web-specific as Django is.

I haven't heard of people using PostGIS to do geoprocessing as much as I see people using QGIS or ArcGIS, so I want to know if it is a comparable alternative.


texts: On the calculus of joint-count test on no-full grid

texts:
I have a data points  of presence-absence records and I want to perform joint-count test as it was suggested in the post https://gis.stackexchange.com/a/388939/178392  by calculating it on a grid (raster). The raster have to contain no-NA values only in the records cells I used for this rasterFromXYZ function from raster package. At this point sometimes I get NA for the p_values of the test. Is it correct what  I did?
cells_data <- cellFromXY(Raster_[[1]], xy =xy ) # Raster_[[1]] is the raster of reference and xy are UTM locations of my presence-absence data
xy_cells <- xyFromCell(Raster_[[1]], cell=cells_data)
data_cells <- as.data.frame(cbind(xy_cells, data$P_A))
colnames(data_cells) <- c('x', 'y', 'P_A')
coordinates(data_cells) <- ~ x + y
proj4string(data_cells)=crs(Raster_[[1]]) #
Raster_PA <- rasterFromXYZ(data_cells, res= res(Raster_[[1]]), crs= crs(Raster_[[1]]))
Raster_PA_ <- extend(Raster_PA, y=extent(Raster_[[1]]))
Raster_PA_ <- mask(Raster_PA_, Raster_[[1]] )
pA <- rasterToPolygons(Raster_PA_, dissolve=FALSE)
nb <- poly2nb(pA, queen = T)
W <- nb2listw(nb, style='B', zero.policy=TRUE)
JCT <- joincount.test(as.factor(pA$layer), W , zero.policy=TRUE)


texts: Select a subset of lidar tiles (las files)

texts:
I have terrestrial lidar scans that were merged and tiled into 100 x 100m tiles, but they are a little too big so I re-tiled to 50x50m. However, a lot of tiles at the edge of the study area are very small and do not contain useful information (e.g., 10x10m) which causes problems in my processing pipeline.
How can I select a subset of tiles from the las catalog picture below?


texts:
I have a vector file with a single, simple polygon (perhaps 10 nodes). I want to add more nodes to this polygon, retaining the original shape.

The obvious question here is 'why would you want to do that?'. It solves an issue further down stream in the process. It's a hack.

Is that possible in FME? Can you read in a vector file, add nodes to it, and write it out again? Which transformer might you use to achieve that?


texts:
I have two layers, one polygon and one of intersecting points. I would like to add a field in my polygon layer that lists the number of points intersecting each polygon.

Original: J'ai deux couches, celle de polygones et celle de points qui s'intercectent. Je veux avoir dans un des champs de ma couche de polygones, le nombre des points respectifs  qui intersectent chaque polygone. Comment faire? Quel outil ou quelle requête utiliser?


texts: getValue() from second order ui.Select in Google Earth Engine

texts:
I would like two make two ui.Select widgets. The items of second ui.Select depends on first ui.Select. Then I would like to get values of both ui.Selects using ui.Button. I have drafted a code, but it does not work as I want. When I click the button, it gets the value of first ui.Select and prints it, but ignores the second ui.Select's value. Is anyone can help me improve the code and give an suggest how can I make it?

https://code.earthengine.google.com/1796c9feaca7638dc195f52850bca64c

var panel = ui.Panel();
panel.style().set({width: '200px', position: 'bottom-right'});
ui.root.add(panel);

var firstname_select = ui.Select({
  items: [{label:'AAA', value:'AAA'}, {label:'BBB', value:'BBB'}, {label:'CCC', value:'CCC'}],
  onChange: function() {

    if (firstname_select.getValue() === 'AAA'){
      var secondname_select = ui.Select({
        items: [
        {label:'A1', value:'A1'}, {label:'A2', value:'A2'}, {label:'A3', value:'A3'}],
        onChange: function() {
        }
      });
    }

    if (firstname_select.getValue() === 'BBB'){
      var secondname_select = ui.Select({
        items: [
        {label:'B1', value:'B1'}, {label:'B2', value:'B2'}, {label:'B3', value:'B3'}],
        onChange: function() {
        }
      });
    }

    if (firstname_select.getValue() === 'CCC'){
      var secondname_select = ui.Select({
        items: [
        {label:'C1', value:'C1'}, {label:'C2', value:'C2'}, {label:'C3', value:'C3'}],
        onChange: function() {
        }
      });
    }
    panel.widgets().set(1, secondname_select);
  }

  });

var getvalue = ui.Button({
  label: 'Get Values',
  onClick: function() {
    print('First: ' + firstname_select.getValue());
    print('Select: ' + secondname_select.getValue());
  }
});

var secondname_select = ui.Select();
panel.widgets().set(1, firstname_select);
panel.widgets().set(2, secondname_select);
panel.widgets().set(3, getvalue);


texts: Selecting several disseminate features with ArcGIS API for JavaScript?

texts:
I have a cadastral map similar to this one from an ArcGIS Javascript API example:

Cadastral map

What I would like to do is selecting several parcels clicking on them and show the information related in a sidebar or floating panel. For example, I would like to select that big house on the right side and then few more spread on the map. Every time I click on one, I would like the application to show the information associated AND not losing that selection when clicking on a new one. In the end, I should have an array of selected features (or similar) that maybe I can re-use later.

I have tried Feature Layer with Selection example, but it does not fully fit my requirements, because depending on the distance between features, I can't select two of them without selecting many more that I do not need.

Is this possible? Can anyone explain the process?


texts: Looping and adding features in memory layer using PyQGIS?

texts:
I would like to know how to add each object obtained as an entity in the same layer; Not by creating a memory layer at  each time in this example.

rb = QgsRubberBand(qgis.utils.iface.mapCanvas())
  rb.setColor(Qt.red)

  geom = QgsGeometry.fromPolyline(p)

for pnt in p:
    rb.addPoint(pnt)

v_layer = QgsVectorLayer("LineString", "cable", "memory")
QgsMapLayerRegistry.instance().addMapLayers([v_layer])
pr = v_layer.dataProvider()
v_layer.startEditing()
v_layer.addAttribute(QgsField("id", QVariant.String))
v_layer.addAttribute(QgsField("type", QVariant.String))

seg = QgsFeature()
geom = rb.asGeometry()
seg.setGeometry(geom)
pr.addFeatures( [seg] )
v_layer.commitChanges()


texts: IRasterLayer background transparency in C++ with ArcObjects 10.x .NET

texts:
I have a georeferenced IRasterLayer which is displaying data from a PPM file with no alpha channel. I am trying to specify a color which should be treated as transparent. In ArcMap it is possible to set a specific color to be transparent, but I have yet to find a working C++ example for ArcObjects 10.x .NET.

I am using the C++ COM interface and have been trying to achieve this by creating a RasterRGBRenderer as per this post. The claimed-working VB.NET sample passes an array of doubles into IRasterStretch.BackgroundValues which is not compatible with the documented 10.1 interface.
Unfortunately documentation in this area is poor, and these samples are very dated. The documentation does not even specify what type of COM VARIANT is needed when VARIANTS are expected.

Please advise how to specify background transparency using C++ in ArcObjects 10.x .NET

Edit:
As a test I have attempted to set the BackgroundValue white to display as red using the following, this did not work (This does work on a test PNG with a white background, but not my PPM file):

esri::IRasterRGBRenderer2Ptr rgbRenderer;
rgbRenderer.CreateInstance(__uuidof(esri::RasterRGBRenderer));
IRasterRendererPtr rasterRenderer(rgbRenderer);
rasterRenderer->Raster = rasterLayer->Raster;
hr = rasterLayer->putref_Renderer(rasterRenderer);

IRgbColorPtr backgroundColor;
backgroundColor.CreateInstance(__uuidof(esri::RgbColor));
backgroundColor->put_RGB(RGB(255, 0, 0)); // Bright red for testing

CComSafeArray<double> backgroundValue;
backgroundValue.Create(3);
backgroundValue.SetAt(0, 255.0);
backgroundValue.SetAt(1, 255.0);
backgroundValue.SetAt(2, 255.0);
CComVariant varBackgroundValue(backgroundValue);

IRasterStretch2Ptr stretch(rasterRenderer);
hr = stretch->put_Background(VARIANT_TRUE);
hr = stretch->putref_BackgroundColor(backgroundColor);
hr = stretch->put_BackgroundValue(varBackgroundValue);
rasterRenderer->Update();

I can confirm that BackgroundValue is being set correctly. Calling get_BackgroundValue always returns a VARIANT array of doubles and these include the correct values if I get after I call put. The renderer is also confirmed to be set correctly as if I disable one of the bands using the renderer pointer that is clearly visible