texts: Precompiled windows binaries for pgRouting latest - v.2.4.1

texts:
I was trying to find precompiled windows 64-bit binaries of the latest pgRouting release Version 2.4.1 but all possible sources like

http://postgis.net/windows_downloads/
http://download.osgeo.org/postgis/windows/pg96/

just list 2.3.2 as the highes available version. Is there any other approach than compiling the latest build on my windows machine?

For linux a precompiled version can be retrieved in form of the package Postgresql-9.6-pgroutingthrough ubuntugis-unstable. I am using the following source:

deb http://ppa.launchpad.net/ubuntugis/ubuntugis-unstable/ubuntu xenial main

Any suggestions?


texts: How to export metadata from a file geodatabase as XML

texts:
I have a file geodatabase with ~20 layers and I would like an XML file for the metadata for each layer that I can send with my data.  I tried the "Export XML Workspace Document" tool but it only gives me one large file instead of individual files.

Is there a way to export individual XML layers of the metadata from a geodatabase?


texts:
I have a layer that was constructed by multiple dissolve, merge, join's etc. by a previous person. So unsurprisingly the table is a complete mess. I only need 2 or 3 out of the 20 something fields that are in the table. However I can't seem to delete them, as the Delete button is grayed out. I found online that there is a workaround using the Delete Fields tool. But that produces:

  ERROR 000840
  The value is not a Raster Layer. ERROR 000840: The value is not a Raster Catalog Layer. ERROR 000840: The value is not a Mosaic Layer.

Which gives dozens of hits here on GIS.stackexchange.com but they vary wildly by topic, but mostly about arcpy or the modelbuilder.

The layer is saved in a layer package, if that might make a difference?


texts:
The below code has written for ASTER radiometric correction in Google Earth Engine but in visualization step, returns this error:

code link:https://code.earthengine.google.com/3ccaaaa0e9ce68bdf80b8056be5d057d
Map.centerObject(table);
Map.addLayer(table)

var ast = ee.ImageCollection("ASTER/AST_L1T_003")
.filterDate('2000-01-01','2020-01-01')
.filterBounds(table)
.filter(ee.Filter.calendarRange(7,9,'month'))
.select('B0[4-9]')
.map(function(img){

  var gain4 = ee.Image(ee.Number(img.select('B04').get('GAIN_COEFFICIENT_B04')));
  var gain5 = ee.Image(ee.Number(img.select('B05').get('GAIN_COEFFICIENT_B05')));
  var gain6 = ee.Image(ee.Number(img.select('B06').get('GAIN_COEFFICIENT_B06')));
  var gain7 = ee.Image(ee.Number(img.select('B07').get('GAIN_COEFFICIENT_B07')));
  var gain8 = ee.Image(ee.Number(img.select('B08').get('GAIN_COEFFICIENT_B08')));
  var gain9 = ee.Image(ee.Number(img.select('B08').get('GAIN_COEFFICIENT_B09')));

  var b4 = img.select('B04').multiply(gain4).clip(table);
  var b5 = img.select('B05').multiply(gain5).clip(table);
  var b6 = img.select('B06').multiply(gain6).clip(table);
  var b7 = img.select('B07').multiply(gain7).clip(table);
  var b8 = img.select('B08').multiply(gain8).clip(table);
  var b9 = img.select('B09').multiply(gain9).clip(table);

  var stack = ee.Image.cat([b4,b5,b6,b7,b8,b9]);
  return stack
}).median();

print(ast);

Map.addLayer(ast)


texts: "Join Attributes By Nearest" in QGIS Graphical Modeler Shifts Attribute Field Order

texts:
I am using the QGIS Graphical Modeler to sample a DEM using two created point features (road & pipeline) then use a spatial join to find the three closest road locations to a pipe location.  The output of the sampled features are as expected, but when I try and use the 'Join By Nearest' tool in the Graphical Modeler to find the closest sampled road locations to sampled pipeline points, it shifts the joined column data of input layer 2 (road) to the right one spot.
However, when I do this same process manually it works out ok.


texts:
I dread installing QGIS because this additionally always installs QT Assistant, QT Designer, QT Linguist. Any reason behind it?
I don't use any plugins but sudo apt install qgis forces to install mentioned QT programs.


texts: Is it possible to load features into ArcSDE without using ArcObjects?

texts:
I have a SDE geodatabase in MS SQL Server, where geometries are stored in MS SQL's native geometry types, rather than the ArcSDE Binary.

It seems like it would be relatively easy to create new features in C#, without the use of ArcObjects.  I could use one of my favorite Micro-ORMS (Dapper or PetaPoco) in conjunction with the Microsoft.SqlServer.Types namespace and load features directly into the database without ArcObjects.

I am a novice ArcSDE user, but am quite aware of the rat's nest of stored procedures and other db logic that comes along with ArcSDE, so I am sure its not as easy as I am thinking.  I am sure I am not the first person who ever thought about this - so the lack of success stories in this area discourage me a bit.

Has anyone ever done this?  Are than any compelling reasons why I shouldn't do this?


texts: Masking raster by using extract by mask tool if feature mask is larger than raster area in ArcGIS Desktop?

texts:
I have a problem to mask the the raster which has values ranges -0.73106217384338 - 0.94593381881714. These values have been derived by using equatin (B3-B4)/(B3+B4) from Landsat band 3 and band 4 spectral reflectance values in arcgis 10. Raster image was of 32 bit tif file formate. I applied the function Extract by mask in arcgis with input this raster and polygon shape file as feature mask. But polygon that i used here cover some more area than the raster. The output raster have some real values like 0.435936 here for the areas where there is no raster beneath the polygon.

How can i get rid of this value using ArcGIS 10?


texts:
I have loaded the following file http://openstreetmapdata.com/data/coastlines to postgreSQL. In visual studio with gdal/ogr and geos enabled:
SELECT geom
FROM waterpoly
WHERE geom IS NOT NULL LIMIT 10;

Stored the result in std::vector* multigeom;

In a for loop

OGRGeometry* geo  =  multigeom->at(i);
printf("geo type: %ld \n",geo->getGeometryType());

I expected it to be 'wkbMultiLineString' (because it says line strings in the above osm website) but the answer was  'wkbPoint'.

Regardless, I have done the following:

OGRMultiLineString* mult = (OGRMultiLineString*) multigeom->at(i);
for( int j = 0; j < mult->getNumGeometries() ; j++ )
{
OGRGeometry* geogeom =  (OGRGeometry*) mult->getGeometryRef(j);
printf("geo type: %ld \n",geogeom->getGeometryType());

}

Answer was 'wkbunknown'

Can anybody give me the reason behind this? or am I doing anything wrong?


texts:
I have point data which has a number of fields/columns, each of which is at a different scale. I'm looking to find a way of constructing a very simple index which would combine this data. Say I have 3 levels of scale (so 3 fields) - small, medium, large. Each of the points either has or hasn't got a text string at each of those levels. I'm looking to assign a number (or not) depending on whether the text string is present or not, and aggregate that data.

Basically I want to see if there are patterns in the presence/absence of text strings at these different scales (possibly re-combining the value given to each of them in more than one way). I'll then figure out a way of showing this map format (Arc). I don't think I've been very clear but please ask if you need clarification.

Here's a sample table with something along the lines of what I was thinking of doing (in grey). I'm not sure this will give me the outcome I need - ideally each column should have a different weighting reflecting relative size. But should give a rough idea of the sort of thing I was thinking of doing.


texts: QGIS crashes when adding ArcGIS FeatureServer Layer

texts:
I want to add an ArcGIS FeatureServer Layer into QGIS. I can add the layer sucessfully but as its being downloaded if I ask QGIS to do anything else it crashes.

First I bring up the Add ArcGisFeatureServer Layer from a Server dialog box. I can connect to the server fine as you can see in the screenshot. The url for the resource is http://maps.six.nsw.gov.au/arcgis/rest/services/public/NSW_Survey_Mark/MapServer. Then I click Add.

At first the dataset appears to be downloading fine.

When I attempt to zoom in on an area it initially looks ok but shortly crashes.

If I choose the debug option I get the following exception.

Am I adding the web service incorrectly? I've tried ticking on the Only request features overlapping the current view extent. No difference. I've added the layer when I was at a zoom level where I knew there are only 3 points. I get the blue loading bar but no features turning up. I assume its downloading them all over the state and not at the extent I'm zoomed to. If I attempt to move around then the same crash occurs.

I'm running QGIS 2.16.3.


texts:
I have used tasklist /svc | find "arc.exe" in Command Prompt, but
in the services.msc list without arc.exe


texts: Previewing raster catalogs raster field in Windows form?

texts:
What I want to achieve:
The attribute table for an ArcGIS Raster Catalog contains a raster field.  When a raster field cell is clicked with the mouse, a preview window displaying the raster associated with the row is displayed.  I want to achieve the same or equivalent functionality, but within a windows form.

What I've tried:
I have managed to display the full attribute table of a RasterCatalog layer within a .Net DataGrid using this code
'Bind a geodatabase table to a .NET control'.   However, whilst the raster column is populated with text 'raster', it is not possible to get a preview of the raster.

One solution could be to add the image for each row to the datagrid as detailed here, but I'd prefer if I could keep the Esri raster field functionality if possible.


texts: ArcEngine 10 Editor:  Split/Rotate/Cut Features toolbar?

texts:
I am currently looking for the Split/Rotate/Cut/Reshape tools in the toolbar control provided by ESRI, however according to ArcGIS Engine commands page there isn't one.

Is it called something else or does it only exist in ArcEditor?


texts: Georeference option not visible in QGIS 3.0 for French version

texts:
I am using QGIS 3.0 Girona Early Adopter Release. I used to use Raster-->Georeference--> ... to georeference images in QGIS 2.18 but now I can only see the interactive georeferencing option in my version. Where did it go?

Could it be because my version is the French release?

When I go to the plugins and search for georef, I don't have the option for georeferencer GDAL only freehand raster georeferencer is available.  I have the experimental plugins option ticked too.


texts: Equivalent functions in ArcMap compared to MapInfo: Objects > Combine

texts: What is the equivalent of MapInfo's Objects > Combine in ArcMap? I have a polyline layer which is divided into 10m subsections, I want to combines these and find parts which are >50m combined. I have done this before using MapInfo by combining the disaggregating the lines (multipart to singlepart), but ideally I'd like to be able to do this in Arc, rather than ping-ponging a file between .shp and .tab format.


texts:
I am using QGIS 3.18 on Windows. Whenever I edit a legend text in Layout view my computer gets very worked up. My CPU temp rises dramatically and the process of typing is really slow. Every single letter that I type takes perhaps 5-25 seconds to type. Instead I am using copy+paste, then I only have to wait once. But a process that could take a minute (to edit 5 rows of legend text) takes much longer time.
Has anyone else encountered this problem?
Edit: Actually the process of editing anything about the legend, for example removing layers that are not to be included in the legend, takes a very long time.


texts: Importing LiDAR (.las) data in PostgreSQL with PDAL?

texts:
I am following this tutorial  to load LiDAR data into PostGIS/PostgreSQL database.

The order is:

Create extensions:

CREATE EXTENSION postgis;
CREATE EXTENSION pointcloud;
CREATE EXTENSION pointcloud_postgis;

Look at the metadata:

pdal info --input 20090429_42122c8225_ld_p23.laz --schema
pdal info --input 20090429_42122c8225_ld_p23.laz --metadata --xml

Run the data loading:

pdal pipeline laz2pg.xml

In this step I have an error message:

  PDAL: laz2pg.xml: JSON pipeline: Unable to parse pipeline:
  * Line 1, Column 1
   Syntax error: value, object or array expected.

Maybe someone knows how to solve this problem?


texts: Ubuntu 20.04 with QGIS 3.14 (pi) SAGA is missing

texts: After updating to QGIS 3.14 in Ubuntu 20.04, SAGA is not available. The Log Panel (View->Panels->Log Messages) said SAGA didn't load. In the Processing Toolbox, the SAGA category is missing.


texts: Does someone have a beginning to end example of how to take data from a SQL Server geography type field and pass it to the Google Maps API to display the shape on a map?  The geography data I'm storing is either a polygon or multipolygon, some of which have holes in them.  I would like to do this using .NET code and JavaScript on the fly with whatever data is stored.


texts: How to use GeoServer SQL views in Openlayers 4

texts:
This is the my script that I am currently using to show my GeoServer layers in OpenLayers:

<!DOCTYPE html>

    <title>WMTS Layer from Capabilities</title>
    <link rel="stylesheet" href="https://openlayers.org/en/v4.6.5/css/ol.css" type="text/css">

The intersect_new layer is a SQL view, which is working and uses default parameters I set in GeoServer. I would like to use the viewparams to change the parameters as in the below link:

http://localhost:8080/geoserver/ui/wms?service=WMS&version=1.1.0&request=GetMap&layers=ui:intersects_new&styles=&bbox=-84702.61914736108,-9272.577651805477,676223.7241900009,1242876.667023777&width=466&height=768&srs=EPSG:27700&format=application/openlayers&viewparams=district:Barnet

I have gone through the OpenLayers documentation and I can't find viewparams as an option for either WMS or WMTS requests. Does anyone know how I can add this as an option to my fetch request?


texts: Upload to Earth Engine from Google Drive on python

texts:
I want to replicate the python keras example on Colab but using Google Drive instead of Google Cloud. I was able to make it work for most of the code instead of the last part consisting of uploading the TFrecord and JSON file to EE, to visualize the results:

!earthengine upload image --asset_id={outputAssetID} {outputImageFile} {jsonFile}

This command line seems to only work with Cloud Storage Buckets and not with Drive. I also tried to manually upload the files to EE but I got the error:

  Manifests for TfRecord ingestion must have exactly one tileset with exactly one source.

I know that it is possible to fix this error using:

!earthengine --no-use_cloud_api upload image --asset_id={outputAssetID} {outputImageFile} {jsonFile}

Is there an alternative to this command line to upload from Drive or a local file from Colab?


texts: Higest-resolution daily temperature products?

texts:
I'm looking to collect daily temperatures for a very small geographic area, and as such am in need of some high-resolution products. I have used MODIS, but the data comes with a 1km resolution.

Are there any alternatives that can achieve the resolution that LANDSAT is able to provide but with more frequency?


texts: Find nearest POI of given GPS coordinates

texts:
I am new to GIS and need to achieve the following:

Find the nearest train station or bus stop of given GPS coordinates in a Java application. Since there are many coordinates to be checked, an offline calculation is preferred.
Visual rendering of maps is NOT required.

After some research, I thought of the following:

use OpenStreetMap, download the maps
extract the relevant nodes using osmosis ("... --tf accept-nodes ...")
import the data into a PostgreSQL/PostGIS database
finally run the queries to select nearby train station

Can you tell me if the outlined approach is possible and useful?
Are there other, better ways to solve the problem?

How can I implement step 4 in a simple yet effective way?
Are there Java libraries which assist with the PostGIS stuff?


texts: r.regression.multi - ERROR: Not enough valid cells available

texts:
I am running a few analyses with r.regression.multi in Grass 7. For a few rasters that I try to correlate, the following error comes up, however: "ERROR: Not enough valid cells available".
I cannot see how the rasters I am trying to correlate in theses cases differ from the rasters where the analysis works. They match in dimension and location.
Any ideas where the problem could lay?

Addition:
Here a screen shot of another situation: the regression is only calculated for a certain area.


texts: Creating equally spaced coordinates on California state map created with basemap?

texts:
I have created a California state map with basemap Python library and this shapefile.

Code is below

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(20,20))
map = Basemap(llcrnrlon=-124.48,llcrnrlat=32.53,urcrnrlon=-114.13,urcrnrlat=42.01,
             resolution='c', projection='lcc', lat_0 =  36.778259, lon_0 = -119.417)
#westlimit=-124.48; southlimit=32.53; eastlimit=-114.13; northlimit=42.01
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='#f2f2f2',lake_color='#46bcec')
map.drawcoastlines()
map.readshapefile('./CA_Counties/CA_Counties_TIGER', 'CA_Counties_TIGER')
plt.show()

Output

Now what I want to do next is to plot equally spaced points only on the California map which would fill the entire California map (no ocean and nearby states). Example would be the grid shown below as all the points in the grid fills the grid(almost).

The way I think it can be done is by considering California as a polygon or multipolygon and generating equally spaced longitude and latitude inside it. I was looking at Generating random coordinates in multipolygon in Python? but it didn't really solve my problem as when I ran the following code to generate points

import fiona
from shapely.geometry import shape
import random
from shapely.geometry import Point
from shapely.geometry import Polygon

def generate_random(number, polygon):
    list_of_points = []
    minx, miny, maxx, maxy = polygon.bounds
    counter = 0
    while counter < number:
        pnt = Point(random.uniform(minx, maxx), random.uniform(miny, maxy))
        if polygon.contains(pnt):
            list_of_points.append(pnt)
            counter += 1
    return list_of_points
all_points=[]
for pol in fiona.open('./CA_Counties/CA_Counties_TIGER.shp'):
    #l.append(pol['geometry'])
    all_points.append(generate_random(50, Polygon(pol['geometry']['coordinates'])))

It gave me the following error

ValueError: A LinearRing must have at least 3 coordinate tuples

Also I am not even sure if the above code worked would it have been given me all the equally spaced points(lat and lon) which would fill the entire California map. Is there other way to do it or someone can help me with the above code. It must generate lat and lon tho. So i can plot them with map.plot() function in basemap. Also would like to get all the points in a datastructure. For example a list of lists or array of arrays or any other if it works well (maybe a dictionary)


texts: Detecting failure of Data Access Insert Cursor?

texts:
I have data that I have now determined is failing because the underlying Oracle table has a unique constraint that fails. When a duplicate row is inserted there are no exceptions thrown - no way currently to know that the operation has failed.

How can I force some kind of feedback? I could do a select / search after doing the insert but that seems clumsy and slow.

    with arcpy.da.InsertCursor(remoteTable, fieldnames) as iCur:
    for val in pCur:
        sr = arcpy.SpatialReference(pval[16])
        pval[16] = arcpy.PointGeometry(arcpy.Point(float(pval[9]), float(pval[10])), sr)
        iCur.insertRow(pval)
    del iCur

If this list has a row that has a duplicate key, it fails and the whole with block fails - however I don't know that it has failed.


texts: MapServer FastCGI connection pool configuration

texts:
I posted this question in the MapServer User mailing list but hoping I can get more eyeballs on it here. I'll update both threads if I make any progress.
The MapServer LAYER documentation says the following on connection pooling in combination with FastCGI:

Additionally, if you have FastCGI enabled, the connection handle will
stay open indefinitely, or according to the options specified in the
FastCGI configuration

I'm using MapServer with PostGIS and FastCGI, and I'm very interested in learning more about this, however I cannot find any information about how connection pooling should be configured in FastCGI.
Can anybody clarify here? How should connection pooling be controlled by FastCGI config?


texts: Loading a .mbtiles map overlays everything with black (as of QGIS 2.16)

texts:
I've been having a frustrating problem with QGIS lately. When I load a raster map with .mbtiles extension in QGIS 2.16, it overlays the whole workspace with black pixels (refer to the attached images).

I tried adding 0 (zero) to "Additional no data value" and it does make the black go away, but at the same time makes the colors in the raster map sort of washed out.

For some reason the new version of QGIS does not properly handle these maps.

Any solutions?


texts: Why is loading qgs/qgz 2 x slower than loading qlr?

texts:
I have a Layer Definition file(.qlr) that takes 50sec to load. When I save it as a project(qgs/qgz) it takes over 2min to load the same data. Why does loading of qgs/qgz take longer than loading qlr?
It seems like the extra time is used after all the layers are loaded(100%), so QGIS seems to be hanging for over 1min.
My data are all based on PostGIS-views, and are both point, lines and polygons. Symbology is saved in the PostGIS-db.


texts: How to merge several google map images into one?

texts:
I want to have a high resolution aerial photo of a catchment area. I took some images from google maps and georeferenzed these and saved them as a layer file. but how can I merge them now that I have one image for the whole catchment area?

The merge tool didn't work. The images are .png s.
I'm sorry, I'm new to ArcGIS and I couldn't find an answer in the web. Also a mosaic layer doesn't seem what I want...

Thanks for your help!


texts:
I am trying to add client-side feature graphics to a feature layer in a map with a roadmap baselayer using the following code, but nothing seems to appear on the map:

function CreateFeatureLayer()  {
    featureLayer=new FeatureLayer({
        source:[],
        objectIdField:"OBJECTID",
        geometryType:"point",
        spatialReference: { wkid: 4326 },
        fields:[ { name:"OBJECTID", type:"oid" }, { name:"name", type:"string" }, { name:"ktypemid", type:"string" }, { name:"kmid", type:"string" }],
        renderer: { type:"simple",  symbol:{ type: "web-style",  styleName: "Esri2DPointSymbolsStyle",  name: "landmark" } },
        popupTemplate: { title: "{Name}"  }
        });

    map.add(featureLayer);

    const data=[{ LATITUDE: 37.6251, LONGITUDE: -119.085, TYPE: "Title",  NAME: "Name" }];

    let graphics=[];
    let i=0,graphic;
    for (i=0;i<data.length;i++) {
        graphic=new Graphic({
            geometry: { type: "point", latitude: data[i].LATITUDE, longitude: data[i].LONGITUDE },
            attributes: data[i]
            });
        graphics.push(graphic);
        }

    app.featureLayer.applyEdits(graphics);
}

Not sure what I'm doing wrong, after modeling an example. I also want to be able to use my own bitmap, rather than an ESRI icon.


texts:
Could someone give me example how to create a route between two points using ol3 and osrm and give me the instructions (ex in 500m right etc).

I found this already of documenation: https://github.com/Project-OSRM/osrm-backend/blob/master/docs/http.md#general-options


texts: PostGIS: How to return polygon nearest the point

texts:
I have two tables in a PostgreSQL / PostGIS database 1  polygon and the other point data
How do I modify the following query so it returns only the smallest polygon nearest the point
The red in the picture

QUERY
select pt.point_id, ply.polygon_id,ply.admin_level
from point_table pt
join polygon_table as ply
on ST_WITHIN(pt.geometry, ply.geometry)
WHERE pt.point_id = 13424428  -- To Test  only

The above returns
point_id    poly_id admin_level
13424428    -1627120    10
13424428    -295351 8
13424428    -180837 6
13424428    -151336 5
13424428    -58447  4
13424428    -62149  2

I need only  the first row for this point in polygon with admin_level 10
EDIT
This Query will form a link table to join the point to it’s nearest parent polygon
It will look something like this, there could be multiply points  in the same polygon
point_id    poly_id admin_level
13424428    -1627120    10
13424429    -1627120    10
13424430    -1627120    10
13424422    -1627121    4
13424423    -1627121    4
13424434    -1627121    4


texts:
I am very new to scripting with ArcMap and I have been doing fairly well with google and youtube... but I've reached a problem.

I am trying to run a script in which country(ies) I specify from a feature are selected and zoomed to. this is so I can preset a template to begin a map composition with automating the first bits.

what works in python is this;

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
arcpy.SelectLayerByAttribute_management("Country_FeatureLayer", "NEW_SELECTION", " NAME = 'Country' ")
df.zoomToSelectedFeatures()

Where I run into a problem is to make it into a tool script and add the "arcpy.GetParameterAsText(0)" within the script, I don't think I am doing it right.  I've tried;

Country = arcpy.GetParameterAsText(0) # earlier in the script
NAME = arcpy.GetParameterAsText(0) # within the selectlayer line...

I'm sure this gets into the lack of education in scripting for me, which is something I am rectifying, anyone got any help for me?


texts: How to replicate the export data right click function for shapefiles using python for batch exports?

texts:
Updated Question! I'm attempting to find/create a python script that replicates the export tool to execute batch exports of featurelayers into specific shapefiles. I want each input featurelayer to be copied and overwritten to a respective output shapefile. I want each output to have a specific name. Using either the FeatureClasstoFeatureClass_conversion() or CopyFeatures_management() functions I can get one featurelayer to output and overwrite a shapefile in a different location. But how to do this for a list of shapefiles?

import arcpy
Mxd = arcpy.mapping.MapDocument("CURRENT")
arcpy.mapping.ListLayers(Mxd)
arcpy.env.overwriteOutput = True

sampleinputList = ["input_a", "input_b", "input_c"] #these are feature layers in the mxd
sampleout_path = "T:/path" #path I'll be batch outputting &overwriting these shapefiles to
samplenames = ["input_a_name", "input_b_name", "input_c_name"]
sampleOutputnamesList = ["T:/path/input_a_name", "T:/path/input_b_name", "T:/path/input_c_name"] # I want each input to correspond to a respective output name
for featureclasses in sampleinputList:
    arcpy.CopyFeatures_management(sampleinputList,sampleOutputnamesList)
#or maybe this might work?
for featureclasses in sampleinputList:
    arcpy.FeatureClasstoFeatureClass_conversion(sampleinputList,sampleout_path,samplenames)

both instances give me a
    Runtime error
    Traceback (most recent call last):
    File "<string>", line 2, in <module>
    File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 2281, in CopyFeatures
    raise e
    RuntimeError: Object: Error in executing tool


texts: Using ArcGIS java script API in local host

texts:
I am new to ArcGIS javascript API and I want to use it in local host mode.
I definitely follow Esri support instruction for using of it. but unfortunately i can not load it.

What do I write instead of [HOSTNAME_AND_PATH_TO_JSAPI] in init.js and dojo.js host when I want to use this library in local host mode?

replacement -> localhost/arcgis_js_api/library/3.16/3.16/


texts: OpenLayers - Access to GeoServer layers with PHP authentication

texts:
I disable all anonymous requests (WMS/WFS/... ) to my GeoServer, but now I can't access to any layer from my OpenLayers app.

I want to know if it is possible authenticate on my GeoServer via PHP or JavaScript before any OpenLayers request.

If yes, what is the best option? Via PHP cURL, JavaScript or proxy (http://docs.geoserver.org/stable/en/user/security/tutorials/httpheaderproxy/index.html)


texts: Rotate, Scale, Translate in OpenLayers 4 and MapServer

texts:
Sometimes, I get spatial data (usually DGNs but I am not sure the format matters) that in the coordinate control sheet a 'Project Datum Factor' is noted. Say: 1.000073496. If I add them as is in ArcMap, the features won't align with the rest of my data. So in order to align them correctly, I use the "Rotate, Scale, Translate" option in the Transformation tab from the layer properties dialog where I enter the quotient of 1 divided by the Project Datum Factor. In this example 0.9999265094094 as shown below:

After doing that, all the layers line up perfectly.
My question is- is there an equivalent workflow I can use in an OpenLayers/MapServer only environment? Not even sure at what stage or level I should try to apply the transformation. The layers are rendered in the OpenLayers map as WMS and/or WFS layers using MapServer. So could MapServer -somehow- apply the transformation and send the tranformed image or vector data to the OL map? Can OL -again somehow- apply this kind of transformation?


texts: Displaying currency format in CARTO Builder scorecard widget

texts:
I am creating a map in CARTO Builder and I would like to have a scorecard widget showing the average income across all points on the map. In the scorecard, I want the average to be displayed with currency format, e.g. $1,000,000.00 instead of 1,000,000.
From what I have gathered, the scorecard widget only takes numeric types as inputs to the operation.
I have tried:

convert the data into string currency format
convert the data into money data type in SQL. Although this shows up in the data window, the money column is not imported with the rest of the data layer in carto Builder, leading me to believe CartoBuilder does not
recognize money type.

Does anyone have a work around to format the scorecard or know what data type will correctly format it?


texts: Draw regularPolygon from left to right in OpenLayers rather than center out

texts:
When I draw a regular polygon in OpenLayers, it draws from the center out (like a circle). I would like it behave more like Handler.Box does (drawing from the UL corner to LR corner.) Is this possible or do I need to figure out how to get Handler.Box to work?

Whenever I use Handler.Box I get this error:

TypeError: feature.geometry.getBounds is not a function
var bounds = feature.geometry.getBounds();

Which is an error in the OpenLayers.js file.

I have all the code for Regular Polygon (all the after insert code and such) so I would prefer to get it to behave how I want. I couldn't find anything in the docs about it, or maybe I didn't understand it.


texts: How to edit file geodatabase in Python using FileGDB Driver?

texts:
I'm trying to use the GDAL/OGR Python bindings to edit a file geodatabase in Python. I followed this good guide which demonstrates how to install QGIS and a special "FileGDB" driver that allows you to edit file geodatabases. This works in QGIS - I can drag in a file geodatabase, specify that I want to open it with the FileGDB driver and then I can edit it.

Now, I want to do this programmatically outside of the QGIS program. I first tried listing all of the drivers in the Python console within QGIS with the following

from osgeo import ogr
import sys

test_gdb = r"C:\Users\5976632\Downloads\geodatabase\GIS\DDS69ff.gdb"
driver = ogr.GetDriverByName("FileGDB")
cnt = ogr.GetDriverCount()
formatsList = []

for i in range(cnt):
    driver = ogr.GetDriver(i)
    driverName = driver.GetName()
    if not driverName in formatsList:
        formatsList.append(driverName)

formatsList.sort()

for i in formatsList:
    print (i)

The results are - (Note that FileGDB driver is listed)

ARCGEN
AVCBin
AVCE00
AeronavFAA
AmigoCloud
BNA
CAD
CSV
CSW
Carto
Cloudant
CouchDB
DB2ODBC
DGN
DXF
EDIGEO
ESRI Shapefile
ElasticSearch
FileGDB
GFT
GML
GMLAS
GPKG
GPSBabel
GPSTrackMaker
GPX
GeoJSON
GeoRSS
Geoconcept
Geomedia
HTF
HTTP
Idrisi
Interlis 1
Interlis 2
JML
JP2OpenJPEG
KML
LIBKML
MSSQLSpatial
MapInfo File
Memory
MySQL
NAS
ODBC
ODS
OGR_GMT
OGR_OGDI
OGR_PDS
OGR_SDTS
OGR_VRT
OSM
OpenAir
OpenFileGDB
PCIDSK
PDF
PGDUMP
PGeo
PLSCENES
PostgreSQL
REC
S57
SEGUKOOA
SEGY
SQLite
SUA
SVG
SXF
Selafin
TIGER
UK .NTF
VDV
VFK
WAsP
WFS
Walk
XLS
XLSX
XPlane
netCDF

Then I wanted to figure out where the version of Python QGIS is using is located so I could use it externally. I did a

print(sys.executable)

and got

C:\OSGeo4W64\bin\qgis-bin.exe

which is not where Python is installed, it's pointing to the qgis.exe. I found what I believe is the Python.exe in the following location

"C:\OSGeo4W64\apps\Python36\python.exe"

I pointed my IDE to use that version of Python and when I list all of the drivers, FileGDB is missing...My question is why when I try to list all of the Drivers outside of QGIS in Python, is FileGDB not listed? Is QGIS setting some special environment when it opens?

If there is an easier way to do this please let me know. I just want to edit a file geodatabase using Python.


texts:
i want to use composer in pyqgis standalone script to create some easy mas using one raster.

my output image is correct i have the image but that code cant print show  the legend in the map,

any idea where i have the wrong in my code for legend?

here the code

from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4.QtXml import *

qgs = QgsApplication([], False)
QgsApplication.setPrefixPath("my/path/", True)
qgs.initQgis()

mapinstance = QgsMapLayerRegistry.instance()
mapinstance.removeAllMapLayers()
raster='raster.tif'
rasterLyr1 = QgsRasterLayer(raster, 'raster')
rasterLyr1.isValid()
rasterLyr1.loadNamedStyle('new.qml')
mapinstance.addMapLayer(rasterLyr1)

#ledend
mapRenderer = QgsMapRenderer()
c = QgsComposition(mapRenderer)
legend = QgsComposerLegend(c)
legend.model().setLayerSet(mapRenderer.layerSet())
c.addItem(legend)

img = QImage(QSize(800,600), QImage.Format_ARGB32_Premultiplied)
color = QColor(255,255,255)
img.fill(color.rgb())
p = QPainter()
p.begin(img)
p.setRenderHint(QPainter.Antialiasing)
render = QgsMapRenderer()
lst = [ rasterLyr1.id() ]
render.setLayerSet(lst)
rect = QgsRectangle(render.fullExtent())
rect.scale(1.1)
render.setExtent(rect)
render.setOutputSize(img.size(), img.logicalDpiX())
render.render(p)
p.end()
img.save('out1.png',"png")


texts:
I have the following Model which is working fine. It's iterating in all features of my cities center layers and adding it as a facily Location. I want to run one at a time. However when it doesn't attend any demand point, because the cutoff I set it breaks my model, instead of running the next city center.

Do you know how to run through all my cities centers without breaking my model if it doesn`t attend any demand point?


texts: Different results for GEOGRAPHY and ST_TRANSFORM query on geo_point field

texts:
I am finding properties that are within 20K distance using ST_DWithin. But when I try
with ::GEOGRAPHY casting the result is different from ST_Transform. Note: SRID of geo_point field is 4326
Indexed Field:
CREATE INDEX properties_geo_point_transform_idx ON properties USING gist (ST_Transform(geo_point, 26986)) WHERE geo_point IS NOT NULL;

Using GEOGRAPHY,
SELECT id,geo_point,ST_AsText(geo_point)
FROM   properties
WHERE  properties.published
  AND  ST_DWithin(geo_point::GEOGRAPHY, ST_SetSRID(ST_MakePoint(8.398522, 47.000084), 4326)::GEOGRAPHY, 20000)
ORDER BY id

Returns 360 rows
Using ST_Transform,
SELECT id,geo_point,ST_AsText(geo_point)
FROM   properties
WHERE  properties.published
  AND  ST_DWithin(ST_Transform(geo_point, 26986), ST_Transform(ST_SetSRID(ST_MakePoint(8.398522, 47.000084), 4326),26986), 20000)
ORDER BY id

Returns 344 rows
When I try to check which points are missing between them using below query,
WITH GEOGRAPHY AS

(
SELECT id,geo_point,ST_AsText(geo_point)
FROM   properties
WHERE  properties.published
  AND  ST_DWithin(geo_point::GEOGRAPHY, ST_SetSRID(ST_MakePoint(8.398522, 47.000084), 4326)::GEOGRAPHY, 20000)
ORDER BY id
),
ST_TRANSFORM AS
(
SELECT id,geo_point,ST_AsText(geo_point)
FROM   properties
WHERE  properties.published
  AND  ST_DWithin(ST_Transform(geo_point, 26986), ST_Transform(ST_SetSRID(ST_MakePoint(8.398522, 47.000084), 4326),26986), 20000)
ORDER BY id
)

SELECT * FROM GEOGRAPHY LEFT OUTER JOIN  ST_TRANSFORM ON GEOGRAPHY.id=ST_TRANSFORM.id
WHERE ST_TRANSFORM.id IS NULL

I found below 16 points are missing on ST_TRANSFORM part.
"POINT(8.40028 46.82024)"
"POINT(8.39473 46.82052)"
"POINT(8.39801 46.8207)"
"POINT(8.42244 46.82132)"
"POINT(8.39879 46.82075)"
"POINT(8.40522 46.82061)"
"POINT(8.42434 46.82157)"
"POINT(8.40572 46.82071)"
"POINT(8.40572 46.82071)"
"POINT(8.20105 46.88213)"
"POINT(8.3985206 46.820629)"
"POINT(8.4231175 46.8212294)"
"POINT(8.3948126 46.820464)"
"POINT(8.3954042 46.8205728)"
"POINT(8.4064802 46.8208733)"
"POINT(8.4043592 46.8202558)"

So my question is why ST_Transform missed some of the data point ?


texts: Network Analyst for Roads, Planes, and Ferries?

texts:
I'm trying to find out if there's a good way to test best travel routes in Alaska that involves land, sea, and air options. I've worked with Network Analyst so I know how to make and analyze a road network dataset, but I'm wondering if there's a way to include air routes and ferry routes as well??

Thinking I could fudge it a little by determining the distance and average speed of the plane and ferry options and draw that as "roads" in the network dataset. Then I'd classify them as different types of transportation.

Any other more elegant ideas?


texts:
I wonder if it is possible using a custom form  for field calculation to be carried out by means of a calculate button.

I have 4 input fields, I enter in two fields manualy a number, after pressing the button the rest of the fields are calculated.

Example:

Pulse: 100-> By user
average:-> calculation of Pulse/3
Steps: 234-> By user
quantity:-> calculation of Average * Steps

Is this possible to make with QT such a form?


texts: Losing CRS by using grass tools in the modeler

texts:
maybe this question has been already answered by somebody. I'm using the graphical modeler alot. A main issue is that some tools (I guess it's the grass tools) are producing results without the correct CRS.

e.g.

Raster (EPSG 3857) -> r.series (average) -> result (EPSG USER 100000).

Is this maybe an issue caused by wrong settings?

Regards.


texts: ArcGIS client user connection hanging when connecting to PostgreSQL geodatabase

texts:
ArcGIS Catalog and Pro are both hanging when I try to connect as a user that is a member of a role. I think it is related to the SDE schema permissions but I'm not sure.

I have an enterprise geodatabase on PostgreSQL 9.6. using ArcGIS 10.6.1.

I have an SDE schema for the geodatabase tables and a data schema, geodata, for feature classes etc. I also have users of the same names.
These users and schemas work fine. I can connect to geodata and create and view/edit FC no problem.

I then created a new role called jpa and a new user called jpa_viewer. jpa_viewer is a member of jpa. I want to use the role to control users. But when I try to connect to the geodatabase using jpa_viewer ArcCatalog or Pro both hang and I have to manually kill them.

Permissions I have set

grant usage on schema geodata to jpa  ;
GRANT select ON ALL TABLES IN SCHEMA geodata TO jpa ;
grant usage on schema sde to jpa;
GRANT select ON ALL TABLES IN SCHEMA sde TO jpa ;

I can connect and query the geodata and sde schema and tables as jpa_viewer no problem using pgAdmin. The problem only occurs using an ArcGIS client.

When I am creating a new connection in ArcCatlog using jpa_viewer I can initially see the list of databases in the instance in the dropdown. It only hangs when I pick one and save.


texts:
I am using a feature layer in JSAPI, where I edit a feature by the standard selection, change attributes, applyEdits process. It seems that there is a bug in JSAPI both 3.12 and 3.14 which makes the apply edits fail on certain features. The callback just errors out without any clue.

The interesting observations:

Failure only happens on certain features and all cases belong to a particular layer (other features in that layer are just fine).
Changing from 3.12 to 3.14 the features that cause the error are changed but did not go away.

Here is a quick snippet of the code:

sq = buildSelectionQuery();
if (sq) {
    all(assetsFeatureLayers.map(function (l) { return l.selectFeatures(sq, esri.layers.FeatureLayer.SELECTION_NEW).promise;})
    ).then(function (fls) {
        console.log('sel res', fls);
        all(fls.map(function (r, i) {
            var fs = r[0]; // the first is an array of selected features
            var l = assetsFeatureLayers[i];
            console.log(fs);
            if (fs.length > 0) {
                console.log("Switching phases to: ", tph);
                fs.forEach(function (f) {
                    f.attributes['phasecode'] = tph;
                });
                console.log("Saving switched phases for layer: ", l.name);
                return l.applyEdits(null, fs, null); //.promise;
            } else {
                return null;
            }
        })).then(function (l) {
            console.log("Phase switching finished successfully!", l);
            clearAllSelections();
        }, function (e) {
            console.log("Error switching phases!", e);
            clearAllSelections();
        });
    });
}


texts: Normalizing and adding kernel density rasters

texts: I have 4 years of road kill data for one road, though each year has a different sample time frame (i.e. May-June vs. April-October) and therefore each year has a different number of surveys and totals of road kill.  I would like to compare kernel density rasters to find the worst kill locations when the 4 years are combined, but I do not want to bias the data toward the year with the greatest number of kills. (I am assuming kill locations in years with fewer samples are representative).  How can I normalize the rasters so I can add them and find the greatest kill hot spots on the road?


texts: Drawing text, callouts, rectangles, etc in QGIS

texts:
In ArcGIS Desktop, you have a specific toolbar for adding text/ drawing rectangles, circles, ... / adding callouts (see screenshot), which I use very regularly.
However, I am switching to QGIS 3.4, since it's free and easy to use. I wonder if there is a similar toolbar in QGIS, like the one in ArcGIS Desktop, because I haven't found it up until now. There are a lot of different plugins that do the job, but it's more complicated than the simple toolbar in ArcGIS Desktop.


texts: Importance of CRS in Getis-Ord Gi* hotspot analysis

texts: I have a question about the importance of the chosen CRS in a Getis-Ord Gi* hotspot analysis - I am running a hotspot analysis (using Zone of Indifference spatial relationship) on a quite restricted dataset around the Galapagos Islands (equator area) and had to change the (projected) CRS from World_Gall to NAD27_UTM_15N. Turns out the hotspot analysis is quite different from the first run to the second. I have read about how the Getis-Ord Gi* works but still am puzzled about how different the results are, especially as the cell area (1 x 1 degree grid) my data is in hasn't changed much. Can anyone explain to me why this is the case?


texts: Rotate OnlineResource in SLD for GeoServer

texts:
How can I rotate an OnlineResource (Raster image fill that can't be auto-generated by QGIS)?
I tried: this answer and also just to write a "Rotate" tag everywhere, it didn't do the job.

Edit: I'm using GeoServer 2.11.2, and this is the sld I'm using:

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor version="1.0.0"
    xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
    xmlns="http://www.opengis.net/sld"
    xmlns:ogc="http://www.opengis.net/ogc"
    xmlns:xlink="http://www.w3.org/1999/xlink"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" >
  <NamedLayer>
    <Name>fill</Name>
    <UserStyle>
      <Title>fill test</Title>
      <FeatureTypeStyle>
        <Rule>
          <Name>test</Name>
          <Title>test</Title>
          <ogc:Filter>
            <ogc:PropertyIsEqualTo>
              <ogc:PropertyName>CODE</ogc:PropertyName>
              <ogc:Literal>1250</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>
          <PolygonSymbolizer>
              <Fill>
               <GraphicFill>
                 <Graphic>
                   <ExternalGraphic>
                    <OnlineResource
                       xlink:type="simple"
                       xlink:href="1250b.png" />
                    <Format>image/png</Format>
                   </ExternalGraphic>
                 <Size>100</Size>
                 <Rotation>45</Rotation>
                 </Graphic>
               </GraphicFill>
             </Fill>
          </PolygonSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>


texts:
I am trying to mosaic NAIP into mosaic dataset and ended up like this with the block.

How do I remove the blocks ?

Does anybody have a suggest for me to fix them ?

I am using ArcGIS 10.1 SP 3 and I also use the Mosaic Dataset.


texts: Seeking 1-band raster colour table workaround in QGIS

texts:
Does anyone know what the workarounds are for the 1-band raster colour table plugin in QGIS?

Currently I find it will run once, then crash on the second iteration. The workaround I use is to delete the two files which the plugin creates in the directory of the raster - *.aux.xml and *.bccPAL1. After this I find I can run the tool again.

Both Wroclaw and trunk.

But is there something I am overlooking, a simpler workaround?


texts: Find CSV points on main roads or close to street intersections, in QGIS

texts:
I work for an advocacy organisation, running a campaign during a by-election, i.e. contained to a single federal electorate.

I've been using QGIS to interrogate how our supporters are distributed across space - mostly looking at how they are clustered, how many are in given postcodes etc - but am still a bit of a newbie when it comes to spatial queries.

We have ~6,000 supporters in the electorate, and I want to get volunteers to ring them and ask them to host yard signs to make our issues really visible. But in order to prioritise them for contact, I would like to rank these supporters' records in terms of how visible their signs would be.

In the screenshot above (which is from our database, not QGIS), each dot represents a supporter. If I'm looking for places to put signs, the dots I've drawn a bodgy red line around are good prospects because they're on main roads. The ones I've drawn a yellow line around are good because they're on an intersection. And the ones I've drawn a [edit] green line around are both: they're at the intersection of one or more main roads.

What I'd really like to do is automatically sort all the records (represented by the dots) sort them according to how significant the street/road is they're on and/or whether they are at the intersection of two streets.

Is there a straightforward way to do this in QGIS?


texts: How to convert raster negetive pixel value to positive

texts: I was working on ground water Level. I used raster math minus function  and the Output raster  grid data contains negetive cell values. ma resulting raster contains values ranges from High : 0,187477 to Low : -1,14927. But I do not want any negetive values.I just need to alter the negetive sign to positive vales.How can I do that using arcgis10.2


texts: Simple user input in QGIS with pyqt plugin

texts:
SplitPolygon plugin allows the user to enter an angle by which the polygon(s) will be split. I would like to modify the plugin and take that value from graphical user input (say, parallel with one side of the polygon). So how does one get simple graphical user input in QGIS? Is there something like Utility in AutoCAD's VBA for points, angles, rectangles?

Just to make it clear, I want the user to select two points and the code should either get the two points or the angle directly.


texts: How to get "p value" for time series data to detect significant change

texts:
I have yearly maximum NDVI data (1982-2008) and extracted all value to Excel. Also calculated slope value for each point during the 27 years. All of you know the slope value represents changes per unit (per year). Now I would like to look for significant changes in NDVI over the 27 years. So I have to get "p value" for each point. How I calculate it?

I tried it in Excel. Firstly I have to estimate t-test, but i don't know which value should be in place of  'mu' for this case! After estimating it I can find p value using TDIST function in MS excel.


texts:
I am trying to add elevation to a 2D point in a QGIS 3 plugin, preferably in a concise and straightforward way. Not finding this easy.

I think part of the problem is that in QGIS 3, fromPoint() has been renamed to fromPointXY(). I am not very clear on how to create 3D point geometry.

I have tried a number of variations without finding one which works and have not been able to find a suitable example online.

Basically what I want to achieve is PointXY + elevation -> Point, which I would expect to be something like this:

z = 1.0 # or whatever

for l in selectedLayers: # for each selected layer
  if l.geometryType() == QgsWkbTypes.PointGeometry: # point layer

    l.startEditing()
    feat = l.getFeatures()

    for f in feat: # each feature

      # Get existing geometry...
      geom = f.geometry()
      pt = geom.asPoint() # returns PointXY

      # Replace with 3D geometry...
      pt3D = QgsPoint(pt.x(), pt.y(), z)
      geom = QgsGeometry.fromPointXY(pt3D) # error, because pt is 3D
      f.setGeometry(geom)

      l.updateFeature(f)

    l.commitChanges()

How do I do this?

For the benefit of others this is now resolved using the answer(s) below.

I am now adding the elevation like this:

f.setGeometry( QgsPoint(pt.x(), pt.y(), z) ) # x, y, z

It turns out that my problem was not so much adding the elevation, as knowing when I had successfully added it.


texts: Compute weight value based on other attributes?

texts:
I have a shapefile with multiple attributes.  I would like to give a "weight" for each value of those attributes in order to get an "ordering" of my features.

My shapefile represents roads, for each road I have informations like, speed limit (30, 50, 70, etc.); condition of the road (good, bad, average, etc.); light (good, average, bad) etc.

My intention is to give to each value a "weight" -> if the speed limit is 70 then value 3, if the speed limit is 30 then value 1 / If the light is good value 1 and if the light is bad value 3 / etc.

At the end I would like a column with the total of "points" each segment of road get.

Can I do this with the field calculator, or do I need SQL query?


texts: How to, if possible, copy data from an ArcGIS REST Service?

texts:
I do not believe this is a copy of Download Dataset From ArcGIS REST Service or Extracting data from ArcGIS REST endpoint? due to the source and type of service.
I would like to copy the features from a government ArcGIS REST Service. It is a MapServer so I am not sure how, or if it is possible. I tried replicating a code I used to copy a FeatureServer but that did not work and instead got an ERROR 999999 which isn't very helpful for diagnosing what is wrong.
The code is as follows:
# Set environment options
arcpy.env.overwriteOutput=True

# Set workspace
ws = r'c:/ws'
wsGDB = os.path.join(ws, 'EnvironmentalContamination.gdb')

# MO Hazardous Waste from DNR
conFields = "https://gis.dnr.mo.gov/arcgis/rest/services/e_start/e_start/MapServer/0"
# Copy fc from rest service
conCopy = "MOContamination"
memCon = arcpy.management.CopyFeatures(conFields, conCopy)


texts:
I am calculating the azimuth between two points whose position is defined in WGS84 coordinates. If I Open Attribute Table -> Open Field Calculator and ask for

degrees(azimuth(make_point(6.020323,47.248690),make_point(6.083703,47.242157)))

I get 95.885 degrees which is inconsistent with the manual calculation atan2(sin(delta_lon)*cos(lat2),cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(delta_lon)), or with the result provided by the French Geographic Institute (geoportail) web site.

However, if following Azimuth calculation, extending lines I
run in the Python Console

p1 = QgsPointXY(6.020323,47.248690)
p2 = QgsPointXY(6.083703,47.242157)
d = QgsDistanceArea()
d.setEllipsoid('WGS84')
d.bearing(p1,p2)*180/3.1415926535

then I get the correct result of 98.584. Could someone enlighten me as to why the first method provides an erroneous result?


texts:
I am a new to maps and related services for viewing it. My goal is to display maps on a website and provide several options to the user for interactivity, like overlaying layers on each other etc. GeoServer seems to the only and most promising option. There are certain things I am not able to understand completely.

How does WMS protocol fit with HTTP? When a browser sends a WMS request to GeoServer, is this request an HTTP request? And when the GeoServer WMS returns a map, is this again an HTTP response?

If both these request and response are HTTP, then why do we need OpenLayers?

I also want to add other common functionalities of a web site, like forms, to my web site . Do I need to use OpenLayers then. Also I am struggling to understand the set up of servers. Will there be a set up of both a web server and GeoServer service running simultaneously on my machine?

How does the flow of information happens then? My web server receives an HTTP request, it requests GeoServer in turn (using OpenLayers) and the GeoServer service returns back to Web Server.

Am I correct?


texts: How to find out the density or percentage of vegetation in each aspect

texts:
I am trying to find out what is controlling vegetation in an area of Saudi Arabia for my masters project and I am currently focusing on the aspect of slopes in the area.
I have a DEM layer and I have used the aspect tool to show me which direction the slopes are facing in. I also have a shapefile feature class of points of vegetation and wish to find out the proportion/percentage or density of this vegetation in each aspect (so how much of the vegetation is only on south facing slopes etc.) but I am struggling to find a way to do this. I have tried to do it with zonal statistics but I don't know how to get it to accept the different aspects as individual zones.

Does anyone have any ideas how I might go about this?


texts: Kernel Density Estimation Weighting Factor

texts:
I don't understand where is this weighting factor coming from, does it need to be in the quakeLyr in the attribute table named with "mag" or?

# Now, we can run the kernel density estimation algorithm by specifying the mag or magnitude field as our weighting factor:

processing.runalg("saga:kerneldensityestimation", quakeLyr, "mag", 10, 0, 0, box, 1, output)

Source of code


texts: Inspect raster value of a pixel to see bits that are on/off in GRASS GIS

texts:
I am working with Landsat-8 C1-Level 2 imagery and I need the information that NASA provides in the pixel_qa file. The information about the content of the pixel is given by turning on a certain bit, plus some eventual other bits to indicate the confidence of the class (e.g. clouds are indicated turning on the 5th bit, bits 6-7 give the confidence value for the cloud, pixels 8-9 the confidence for the cirrus). This means that there are several combinations of values that can indicate the same class.
I know r.mapcalc has some bit-wise operators, but I haven't found what I'm looking for: is there a way to check if a specific bit (e.g. the 6th out of 16) is on or off?
Is there any more practical way to retrieve the information I need?


texts: Have NDVI layer and want to Calculate area(pixel area-reduce region) .in Python API for Earth Engine

texts:
I had already asked this question for Earth Engine Javascript.The link is given below
I have NDVI layer in Earth Engine. I want to assign and calculate area of good and bad ndvi area
 I wamted to do the same thing in Python API for earth engine. But it is showing the following error-
EEException: Invalid argument for ee.Reducer(): ({'geometry': [[74.022029, 20.103245], [74.022029, 20.116564], [74.038048, 20.116564], [74.038048, 20.103245]], 'reducer': , 'scale': 10},).  Must be a ComputedObject.
I am also sharing a snippet of my code in Python-

    ndvi = clip.normalizedDifference(['B8', 'B4'])
    #NDVI=ndvi.multiply(100).uint8()
    bad=ndvi.gte(0)and(ndvi.lte(0.2)).rename('bad')
    avg=ndvi.gt(0.2) and(ndvi.lte(0.35)).rename('avg')
    good=ndvi.gt(0.35)and(ndvi.lte(1)).rename('good')
    ndvi.addBands(['bad', 'avg', 'good'])
    print (ndvi)
    areas = ndvi.select(['bad', 'avg','good']).multiply(ee.Image.pixelArea()).reduceRegion({'reducer':ee.Reducer
    .sum(),'geometry':geometry1,'scale':10});
    print (areas)


texts: Format of points in Polyline shapes in a shapefile?

texts:
I am have been asked to write some code that looks at a SHP file that contains a load of Polyline shapes that are used to represent parking areas in a British council's area. I've created a C# program to read the data, but I have noted a number of things that don't seem to make sense to me so I wondered if anyone would be able to shine some light on it (I am not a GIS guy)...

Looking at the points on a Polyline shape, it seems as though the points start with extreme numbers for the X and Y coordinates ( -1.4600101262233409E+239 and -8.5448987622302689E-70 for instance) followed by more sensible figures like 185747.93946064217 which seems to bear a relationship to the GPS position (-0.18ish).

Am I right in thinking that this positional data is actually subject to some sort of base box, and if so, where are it's settings stored? Also, how do we set the type of coordinates used?


texts: Can we insert Legend with only features that are in a shapefile?

texts:
I want to insert legend which have 40 labels data, since, the imported symbology have 90 labels, all the labels are getting inserted in the legend.

I do not want to remove them manually from the legend so is there an automated way?


texts:
I am using QspatiaLite in QGIS 2.14.
I am having problem updating a column of my table and I suspect that the types (integer, real, text etc.) are not consistent.

How can I know which are the types of the fields of my table imported from QGIS into QspatiaLite?


texts:
I'm trying to analyze nighttime lights of a specific area with VNP46A1 data. I got my data from this website, but I could find only few tutorial about this product, and one of them was VIIRS_SR_Tutorial. When I listed VNP46A1 Metadata, I found my first challenge.

from h5py import File
data_path = 'path-to-my-h5-file'
file         = File(data_path)
fileMetadata = file['HDFEOS INFORMATION/StructMetadata.0'][()].split()
fileMetadata = [m.decode('utf-8') for m in fileMetadata]
print(fileMetadata[:11])

>>>'GROUP=SwathStructure',
   'END_GROUP=SwathStructure',
   'GROUP=GridStructure',
   'GROUP=GRID_1',
   'GridName="VNP_Grid_DNB"',
   'XDim=2400',
   'YDim=2400',
   'UpperLeftPointMtrs=(20000000.000000,50000000.000000)',
   'LowerRightMtrs=(30000000.000000,40000000.000000)',
   'Projection=HE5_GCTP_GEO',
   'GridOrigin=HE5_HDFE_GD_UL'

In this page, It says "VIIRS/NPP Daily Gridded Day Night Band 500m Linear Lat Lon Grid Night" but in Production Information section it says the data has 750m spatial resolution. However, when I checked this value by dividing sides by pixel counts I found this

(30000000.000000 - 20000000.000000) / 2400
>>> 4166.666666666667

(50000000.000000 - 40000000.000000) / 2400
>>> 4166.666666666667

what is the real spatial resolution of my data ?

Also another challenge for me is was projection types. In the tutorial, the data was clearly projected Sinusoidal coordinate system,that can be seen in metadata of VNP09GA

Projection=HE5_GCTP_SNSOID

and the author defined it like this

prj = 'PROJCS["unnamed",\
GEOGCS["Unknown datum based upon the custom spheroid", \
DATUM["Not specified (based on custom spheroid)", \
SPHEROID["Custom spheroid",6371007.181,0]], \
PRIMEM["Greenwich",0],\
UNIT["degree",0.0174532925199433]],\
PROJECTION["Sinusoidal"], \
PARAMETER["longitude_of_center",0], \
PARAMETER["false_easting",0], \
PARAMETER["false_northing",0], \
UNIT["Meter",1]]'

But in my data, it says the projection of VNP46A1 is

Projection=HE5_GCTP_GEO

In this Manual, section 8.3.1 (GCTP Projection Codes), I found out what HE5_GCTP_GEO stands for, which is Geographic. But still don't know how to define it.

If the the projection, used to define my satellite data, was Geographic, are my projection coordinates and latitude-longitude values are the same? In what coordinate system my corner points are defined?

UpperLeftPointMtrs=(20000000.000000,50000000.000000)
LowerRightMtrs=(30000000.000000,40000000.000000)

I need to georeference my data but with all these unknowns I cannot do that


texts: Get GeoJSON or simple WKT from a GPFeatureRecordSetLayer ArcGIS js API

texts:
Is it possible to get the GeoJson representation of a GPFeatureRecordSetLayer that is returned from GPServer task?

I am able to draw the feature on the map (using the ArcGIS js API) with code like this:

const polygon_feature = data.value.features[0];
polygon_feature.symbol = fillSymbol;
graphicsLayer.add(polygon_feature);

where data is a returned promise from the .getResultData() method implemented in the ArcGIS API.

Is there any method available to get the geometry from this object so I can pass it to my backend and save it to our database as GeoJSON or wkt?


texts: How to georeference multispectral raster without loss of data

texts:
I am trying to align two rasters made of data captured during same flight, but with separate cameras (drone - onboard, and Micasense RedEdge) and two separate gps devices (drone and another one feeding data to RedEdge).

Since both cameras are taking pictures independently and both gps chips feed data with their own error the result is a mismatch of overlaid rasters - clearly visible by using vector lines indicating rows.

My goal is to overlay the rasters in a way that vector layer indicating rows would have same position on both maps. I do not need absolute accuracy therefore in my trials I used data from RGB raster to georeference index map. But after doing so - I loose all values in the index map receiving only black raster with range of values "nan - 0"

How to set up the georeferencer (transformation settings etc. ) to avoid this such situation?


texts:
I am trying to create a series of layers using a county map of the US. In order to save load time I am only loading the county polygon geometries once, then using applyEdits (https://developers.arcgis.com/javascript/3/jsapi/featurelayer.html#applyedits) to add the county geometries to the layers, like a template. I have tried this in a lot of different ways here is what I currently have:

self.initializePolygonLayers = function(){
    var layers = self.layers;

    $.getJSON("api/layerData/templates/counties", {}, function (data) {
        self.initializeSinglePolygonLayer(self.layers()[self.layerIndex], data);
    });
};

self.initializeSinglePolygonLayer = function(layer, data) {

    if(layer.baseLayerName != "County")
    {
        self.layerIndex++;
        if(self.layerIndex < self.layers().length) {
            self.initializeSinglePolygonLayer(self.layers()[self.layerIndex], data);
        }
    } else {
        var counties = new Array();
        console.log("defining layer: " + layer.name);
        window[layer.name] = defineFeatureLayer("ObjectID", 2, false, ["${NAME}", "${*}"]);
        window[layer.name].setRenderer(createRampRenderer(layer.id));
        window[layer.name].setVisibility(false);
        map.addLayer(window[layer.name]);

        for(var index in data) {
            var county_obj = data[index];
            var geometry = wktToEsriJson(county_obj.geometry);
            var attributes = { "NAME" : county_obj.NAME, "ObjectID" : index, "PostalID" : county_obj.Postal_CountyRegion_ID
, "mappedField" : 0.0};

            if(geometry !== false) {
              var graphic = new esri.Graphic(geometry);
              graphic.setAttributes(attributes);
              counties[counties.length] = graphic;
            } else {
              console.log("Error Loading WKT: " + county_obj.geometry);
            }

        }

        console.log("applying edits:" + layer.name + " counties count: " + counties.length);

        window[layer.name].applyEdits(counties, null, null, function(data){
            console.log("success");
            console.log("success for layer " + layer.name + " graphics count:" + window[self.layers()[self.layerIndex].name].graphics.length);
        }, function(data){
            console.log("fail");
            console.log("fail for layer " + layer.name + " graphics count:" + window[self.layers()[self.layerIndex].name].graphics.length);
        }).then(function() {
            console.log("then for layer " + layer.name + " graphics count:" + window[self.layers()[self.layerIndex].name].graphics.length);
            self.layerIndex++;
            if(self.layerIndex < self.layers().length) {
                self.initializeSinglePolygonLayer(self.layers()[self.layerIndex], data);
            }
        });
    }

};

Here is what happens, applyEdits succeeds for all the layer (I'm trying it for two layers) the first layer renders fine and the second one is blank. Looking into it the the first layer gets all the graphics and the second one, despite applyEdits succeeding, has no features. What's going on here?


texts:
Can we find the event when basemapGallery is completely loaded in arc GIS JS API 4.4 ?
Because when try to find the items in basemapGallery view model, it is taking sometime to get it loaded, currently am achieving it through  setTimeout method, kindly suggest if you have any good idea


texts: I created a kml using the kmllayer. When I set the map with a basemap and use the addLayer() method the kml is not shown. I can see the kml layer without a basemap. Can someone post sample code on using basemap with kml? am using ArcGISTiledMapServiceLayer as basemap


texts: Method asPoint().x() doesn't work in QGIS 3.4.3 but did in 3.4.1

texts:
I have built a plugin which does essentially this
equal sized polygons along line
Main parts of the code are from user lejedi76

I have used this plugin for months now and it worked as expected. Turning from QGIS 3.4.1 to 3.4.1 I get the following error:

  Null geometry cannot be converted to a point

I know what the error means, but I do not know why it even shows up. I also could not find the critical changes in the QGIS version that caused this problem. It seems I am overlooking something very simple, but what?
The code provided works without any problems in the QGIS 3.4.1 Python console:

width = 142.365
height = 79.842
overlap = 0.2
gridPath = r'C://qgistesting/testgrid.shp'

#to test take a line-layer with one feature only
layer = iface.activeLayer()
layer.selectAll()
for feature in layer.selectedFeatures():
    geom = feature.geometry()

    pages = QgsVectorLayer("Polygon?crs=epsg:"+str(25833), layer.name(), "memory")

    fid = QgsField("fid", QVariant.Int, "int")
    attributes = [fid]
    pages.startEditing()
    pagesProvider = pages.dataProvider()
    pagesProvider.addAttributes(attributes)
    curs = 0
    numpages = geom.length()/(width)
    step = 1.0/numpages
    stepnudge = (1.0-overlap) * step
    print(3)
    pageFeatures = []
    r = 1
    while curs <= 1:
        startpoint =  geom.interpolate(curs*geom.length())
        endpoint = geom.interpolate((curs+step)*geom.length())
        ##### Null geometry cannot be converted to a point ######
        x_start = startpoint.asPoint().x()
        y_start = startpoint.asPoint().y()
        x_end = endpoint.asPoint().x()
        y_end = endpoint.asPoint().y()
        currline = QgsGeometry().fromWkt('LINESTRING({} {}, {} {})'.format(x_start, y_start, x_end, y_end))
        currpoly = QgsGeometry().fromWkt(
            'POLYGON((0 0, 0 {height},{width} {height}, {width} 0, 0 0))'.format(height=height, width=width))
        currpoly.translate(0,-height/2)
        azimuth = startpoint.asPoint().azimuth(endpoint.asPoint())

        currpoly.translate(x_start, y_start)
        currpoly.asPolygon()
        page = currpoly
        curs = curs + stepnudge
        feat = QgsFeature()
        feat.setAttributes([r])
        feat.setGeometry(page)
        pageFeatures.append(feat)
        r = r + 1

    pagesProvider.addFeatures(pageFeatures)
    pages.commitChanges()

    QgsVectorFileWriter.writeAsVectorFormat(pages, gridPath, "utf-8", pages.crs(), "ESRI Shapefile")
    layer_loaded = iface.addVectorLayer(gridPath, "Grid", "ogr")

And here comes the lazy answer for my problem after the problem was identified by ndawson (lazy answer kills the last polygon, there will be a better answer)

while curs <= 1:
    curs = curs + stepnudge
    if (curs+step) <=1:
        startpoint =  geom.interpolate(curs*geom.length())
        endpoint = geom.interpolate((curs+step)*geom.length())
        x_start = startpoint.asPoint().x()
        y_start = startpoint.asPoint().y()
        x_end = endpoint.asPoint().x()
        y_end = endpoint.asPoint().y()


texts: How to modify a cell value in a raster under certain conditions (ArcGIS)?

texts: I have a DEM raster with some negative values. After overlying other shapefiles, I realize that the negative values are for bodies of water (mainly lakes and bays). I want to replace the negative values with the positive values elevation values (lakes). Any suggestion to do this in python? (Con tool will not work: lakes are at different altitude levels).


texts: ArcMap 10.3.1 Crashes when trying to change the version

texts:
I'm using ArcMap 10.3.1 when I have a data connection to another user (For example change from Tom's version to Jim's version of the SDE) or default and attempt to change the version to my own ArcMap crashes immediately. When selecting the debug report it also crashes. Any ideas?

Windows 7,
ArcMap 10.3.1

Tried the following:

Restart computer
Repair ArcMap install
Delete Normal.mxt file
Reset ESRI user profile

[EDIT - 8/3/2016]
Issue is back previous fix does not seem to be working. Noticed that if log into another computer with my version I can change the version without issue. So points to something wrong with either this install of ArcMap or the computer, but not on the versioning side of things.


texts:
I would like to define my own gdistance::transition function. However, there are two issues that I aren't clear to me yet.

First, I do not yet understand what the input vector to this transition() function looks like. Transition layers based on elevation often use x[2] – x[1]. Others use mean(x), max(x) etc. I assume that the elements of vector x are cell values from the raster layer conditional on the directions option. But which cell is at which position? And what would that function look like if the cost of travelling through a cell only depends on the cell’s value itself, independent of neighbouring cells. Would the function just be x?

Second, how can I feed multiple raster layers into the transition function? Assume that I would like to compute the transition costs between cells based on elevation and temperature. Is there a way to account for the two variables in a single transition function? The documentation indicates that I can use multiple raster bricks but then I don't know how to access the specific elements in my transition function.

The package documentation is illustrative in many regards but unfortunately remains scarce concerning the definition of the actual transition functions. I'm really stuck here.


texts: what I am trying to do is to generate a heat map depicting the spatial distribution of certain species found in a lake, in association with their sex. For example, species_a+male or species_b+female or species_a+female and so on.Also, I will be assigning unique arithmetic values to every combination, e.g. species_a+male=1, species_b+male=2 and so on, in order to use this non-spatial data. So, the question is how can I do it? I already experimented with heatmaps or the simpler Topo to Raster method but nothing seems to do the trick.


texts: How can I make a vertical symbol for a point which stands upright on the centroid?

texts:
Context: I need to display the location of wind turbines using a symbol from a True-type-font(ttf). The GIS software doesn´t allow me to define offsets, so I need to define a ttf which displays a turbine symbol which stands on the coordinate of the point. This is however more difficult than expected, because the font character is always centered on the point, however I design it. This means that the point location is always half way up the wind turbine and not at it´s foot.

Has anyone designed a ttf and found a solution to this problem?


texts: Symbology forMultiple Attributes in ArcMap?

texts:
I need some help to map landslide proximity to urban areas and their size on ArcMap 10.5.1. I am unsuccessfully trying to use symbology by multiple attributes. The data is point data but within the attribute table I have information about landslide proximity to urban areas and their size.

I would like a different colour for the points according to their proximity to urban areas - eg for those less that 1km distance, those more than 1km but less than 2km etc. I would also like the size to be linked to the size of the landslide.

Could someone advise me how to create the ranges for colour section - I can only seem to colour every point arbitrarily.

Unlike for symbology by quantities there seems to be no way to put in a range for the value (see fig). How does one do this?

I wish for a range such as this to be adopted but then the symbols to also be graduated in size according to a field called shape area.

I am unable to change the size and keep different colours for the landslides according to their proximity (as I did above in fig2).


texts: Filtering layers in web map legend based on viewing extent using ArcGIS Online

texts:
I am creating a web map for my organization which has several natural preserve sites in a given state. This is intended as a simple introduction for the public where they can view basic details of all the reserves from one public map i.e. boundary, walking trails, kayak trails as some are coastal.
I'm using visibility ranges such that at the statewide level users only see locations points for our sites and offices, and when they zoom into a specific preserve, trail layers and other specifics appear. Problem is that the legend will show all the layers within the visibility range, even if they're not in the map extent. Ex. I zoom to Site A, there's Site A layer 1, Site A layer 2, but also layers for Sites B, C, and so on. Of course for trails I could merge them into one statewide trail file, but we have various site-specific projects or features that aren't uniform for all sites. Ex. Site C has prohibited federal property but most others don't.
Is there a way to filter content on a web map legend based on what's viewable in the user's extent?
Preferably with a web map since we embed those in our website and want it to be as simple as possible for public viewing.
I have tried using a feature layer view definition, but that defines what features are viewable from that feature layer, not what extent the feature layer can be viewed in by a user.


texts: Error on showing .js as spatial object with JSON.parse()

texts:
I have a Leaflet map exported by qgis2web plugin. The geodata is stored as (Geo-)JSON-objects inside a variable saved as .js-file.
I have an attribute field ("bildpfad") with comma-separated pathes to images I want to show in popups. For showing the images, I need access to every single path. That works by converting the path-field to an array and going through every single element with array.forEach().
Actually I am using the JSON.parse() function. I have other workspaces where this way works great, but not in that specific one.

I got following error message in the browser console:

SyntaxError: JSON.parse: unexpected character at line 1 column 2 of the JSON data

I've already validated the JSON-object without any errors. For testing purposes I reduced the hole dataset to one feature, but with the same error-message. Here  the JSON-data I'm using:

var json_spielgeraete_freizeit_2 = {
    "type": "FeatureCollection",
    "name": "spielgeraete_freizeit_2",
    "crs": {
        "type": "name",
        "properties": {
            "name": "urn:ogc:def:crs:OGC:1.3:CRS84"
        }
    },
    "features": [{
        "type": "FeatureCollection",
        "name": "spielgeraete_freizeit_2",
        "crs": {
            "type": "name",
            "properties": {
                "name": "urn:ogc:def:crs:OGC:1.3:CRS84"
            }
        },
    "features": [{
        "type": "Feature",
        "properties": {
            "bildpfad": "\"U:\\ER-GIS-Verbunddokumente-Baumgrün\\spielgeraete\\1010443.JPG\"",
            "kategorie": "Wippgerät"
        },
    "geometry": {
        "type": "MultiPoint",
        "coordinates": [
            [10.9846280196624, 49.572284008148969]
        ]
    }
}]
}]
}

And here the corresponding JavaScript part:

var string = feature.properties['bildpfad'] !== nullAutolinker.link(String(feature.properties['bildpfad'])) : '';
var path_array = JSON.parse("["+string+"]");
var popupContent =
        '
\

\

\
            test	\

\

\

\
                test	\

\''
                path_array.forEach(path_array => popupContent += '
<img src="'+path_array+'" width="250px"
')
                popupContent += '

';

        layer.bindPopup(popupContent, {maxHeight: 400});
    }

I don't think the umlauts ä,ö,ü causing the problem, removing them doesn't make any difference. The same by adding/removing / or \ in front of unsupported characters like /. I'm new in JavaScript, so maybe I'm misunderstanding some syntax rules in JSON-files (e.g. keywords or quotations marks) or somewhere else.


texts: How to move position of identify features window in QGIS 2.18.10

texts:
When I select the "Identify Features" tool in the attribute toolbar, and then right click on a WMS layer, then the top of the appearing window is outside my screen and I cannot move it around.
How can I get this window to appear in a way, that it can be moved around?


texts: Showing items from one layer based on values in another layer

texts:
I have a few points layers each showing the locations of particular facilities. I have two layers with polygons buffering from centre and from boundary of a site. I would like the correct buffer to appear should a certain set of points be shown on a map i.e if nursery points are shown on map then only the 1Km buffers will show, if I replace the nursery layer with schools layer, the 1Km buffers are removed and replaced by the 2Km buffers.

I think this would be controlled by 'data derived' expressions but I don't know where to start!


texts: os.system /subprocess.call error handling when error in command line function?

texts:
In my Python script I call a GDAL utility using:

subprocess.call('gdal_rasterize -a ICE_TYPE -where \"ICE_TYPE=\'Open Water\'\" -b 1 -burn 5 -l ' + shapefileshortname +' ' +  shapefile + ' ' + outraster)

or

os.system('gdal_rasterize -a ICE_TYPE -where \"ICE_TYPE=\'Open Water\'\" -b 1 -burn 5 -l ' + shapefileshortname +' ' +  shapefile + ' ' + outraster)

Now gdal_rasterize gives an error for a corrupt shapefile, but using "try-except" on os.system only works for errors related to os.system itself.

Is there within the Python script a way to handle errors by a command line function called via os.system or subprocess.call ?


texts: UpdateCursor can't acquire lock, but can acquire lock in TestSchemaLock

texts:
import arcpy
import os.path
arcpy.env.workspace = "C:/Users/melis/Documents/HenryWork.gdb"
arcpy.env.overwriteOutput = True
buffer = "LS_buffer"
landcov = "ForestLC_clip"

#create search cursor to access table. may need to change field name with other data.
buffCursor = arcpy.da.SearchCursor(buffer, ["buff_id"])

for row in buffCursor:
    #select one buffer
    query = """ "buff_id" = %s """ %row[0]
    arcpy.SelectLayerByAttribute_management(buffer, "NEW_SELECTION", query)

    #define clip out path
    clip_out = "LS" + str(row[0])

    #clip w one buffer
    arcpy.Clip_analysis(landcov, buffer, clip_out)

    #add fields to clip_out
    arcpy.AddField_management(clip_out, "clip_id", "TEXT", 20)
    arcpy.AddField_management(clip_out, "area_ha", "FLOAT")

    #create updatecursor for new clip_out
    clipCursor = arcpy.da.UpdateCursor(clip_out, ["clip_id", "area_ha"])

    for polys in clipCursor:

        #calc geom for area ha
        exp = "!SHAPE.AREA@HECTARES!"
        arcpy.CalculateField_management(clip_out, "area_ha", exp, "PYTHON_9.3")

        #field calc for buffer id
        polys[0] = row[0]

del clipCursor

The code works all the way down to the first calculate field, and then I get the error: "Cannot acquire a lock".
However, when I test whether I can acquire a lock on that feature class (the new clip_out ("LS1"), I get a "True".

arcpy.TestSchemaLock("LS1")
True

Any idea why the cursor code can't acquire a lock, but when outside the code, I can acquire the lock?


texts: How do you label only the odd or even integers of a field using python in the label expression in arcmap?

texts:
I'm new to python and would like to use the label expression box with python to display only the odd or even integers from the field, [Frame_Number].

I was thinking something like

def FindLabel ( [Frame_Number] ):
  for num in [Frame_Number]:
     return num % 2
  return [Frame_Number]

But I know that is incomplete.


texts: Auto delete 5th record while inserting 6th record of polygon to SpatiaLite table

texts:
I am having problem while inserting 6th record to SpatiaLite table it auto deletes 5th record, and after that it is ok.

I am using SpatiaLite database to insert some polygon using Google Maps API in an Android application.

My insertion code is like:

      String polyPoints = "POLYGON((";
        if (listofpoints.size() != 0) {
            for (int i = 0; i < listofpoints.size(); i++) {
                polyPoints += listofpoints.get(i).longitude + " " + listofpoints.get(i).latitude;
                if (i < listofpoints.size() - 1) {
                    polyPoints += ", ";
                }
            }
            polyPoints += ")) ";
        }

        try {

            Stmt stmt01 = null;
            stmt01 = database.prepare("INSERT INTO " + table + " (id, username, cutarea, createdate, rosen, range, treecat, fields, GEOMETRY) VALUES (NULL, ?, ?, ?, ?, ?, ?, ?,  ST_GeomFromText('" + polyPoints + "', 4326));");
//                stmt01 = database.prepare("INSERT INTO " + table + " (id, username, cutarea, createdate, rosen, range, treecat, fields, GEOMETRY) VALUES (NULL, ?, ?, ?, ?, ?, ?, ?,   MakePoint(" + listofpoints.get(i).longitude + ", " + listofpoints.get(i).longitude + ", 4326));");

            stmt01.bind(1, username);
            stmt01.bind(2, cutarea);
            stmt01.bind(3, createdate);
            stmt01.bind(4, rosen);
            stmt01.bind(5, range);
            stmt01.bind(6, treecat);
            stmt01.bind(7, fields);

            stmt01.step();
            stmt01.close();
        } catch (Exception e) {
            e.printStackTrace();
            return false;
        }

but while I try to insert polygons in the loop it is ok.

What should I try next?


texts: I'm trying to change the symbol from 'marker' to 'fill' in layer properties.  Is it possible to change that so I can fill the entire property instead of having points for each property?


texts: ArcMap - Blocky/Shattered Glass DEM Fix?

texts: I am currently working with some data obtained from the Spanish Government to produce Digital Elevation Models for a site I am studying. I want to create a model that is essentially showing the changes in elevation smoothly, but when using the 5m DEM's provided, I am faced with a "terracing" (see image) effect upon the maps. Is there a way to remove this such that the data are more continuous?


texts:
I have a line and some points along with line, but they are not equally distanced from line. I want to align them at equal distance along the line,
let's say, every point should be at a distance of 10 meter from the line.

I want to make an Python script for this, but don't know where to start.


texts: OpenLayers plugin, project CRS and grids in Composer

texts:
Apologies in advance for the potential for this to be a duplicate question, I've had a good look, and I can't see a solution that I understand well enough to help me.

I have my project CRS set to a local UTM projection, however when I add a map from OpenLayers the project CRS changes to EPSG:3857.  When I add a grid in Composer, I need to display the longs and lats of the local UTM, not the OpenLayers coordinate reference system.

I understand that I need to have the project CRS the same as the OpenLayers CRS, but is there a way I can:

save and reproject the section of the OpenLayers raster into my local UTM, so it still lines up with my vector layers or
display the appropriate coordinates on the grid, despite the project CRS being in EPSG:3857.

This is an issue that I have with many of my maps, where I use aerial imagery for my base layers, but need to display local coordinates in my final map.

Will there be an easy solution to this in QGIS 2.0?


texts:
Currently I have a problem while developing an app - I use OpenStreetMap and draw a Circle and pass it to a Python back-end, as central point and point on the perimeter on the Circle. Next I am describing that circle with 20-sided Polygon, and return both to the front and I see, that actually is no more circle-like, but it looks like an Ellipse.

More to that - I have written an algorithm which defines if my route crosses that zone - and it works on the data from Polygon, but not the Circle from the front (look picture)

The relevant code is:
Getting Radius from central Circle point and point on the perimeter:

def edit_radius(point, radius):
    xc, yc = point[0], point[1]
    rad_x, rad_y = radius[0], radius[1]
    return math.sqrt((xc - rad_x) ** 2 + (yc - rad_y) ** 2)

The similar result I can reach with mpu.haversine_distance(point, radius)/100 - which returns Haversine distance between that points (and it differs with Euclidean distance).

Drawing polygon with 20 sides, and radius passed from previous function:

def polygon(sides, radius=1, rotation=0, translation=None):
    one_segment = math.pi * 2 / sides

    points = [
        (math.sin(one_segment * i + rotation) * radius,
         math.cos(one_segment * i + rotation) * radius)
        for i in range(sides)]

    if translation:
        points = [[sum(pair) for pair in zip(point, translation)]
                  for point in points]

    return points

I can't understand what I am missing. Should I convert to another coordinate system? But as I have read OpenStreetMap (and Leaflet) are using "flat" EPSG:3857 system.

Can somebody explain what is going on?


texts: Is it possible to run Points2one in batch mode?

texts:
Is it possible to run Points2one in batch mode (using Qgis 2.8.1 on Unbuntu 12.04)?

Basically I have 85 separate point layers and I want to create polygons using Points2one for each layer.


texts: Creating points along line according to attribute data using QGIS?

texts:
I would like to create equidistant points along lines in QGIS, but each of my lines have an attribute of how many points should be on them in a column called "piece" and these values are different in almost every case.

My first thought was to divide lines to equal pieces according to the values in the attribute table and then put points to the centroids of the lines but I could not do it.

This question is very similar to what I would like to do, but it is not for QGIS:

How to create the same number of points along multiple polylines?

Most similar questions that can be found here at stackexchange in this topic are about a specific distance between points, a number that can be easily typed into a field (like in the case of Create points along lines or Convert lines to points algorithms) and not about specific data from an attribute table, which I've been unable to find how to use for such purposes.

It would be fantastic if there were no points right at the start and end point of the line but I would be just over the moon if there were as many points on a line as in its attribute table.


texts: How to add a bounding box filter to this leaflet WFS request?

texts:
I use this function to add an interactive layer to my leaflet application. It works great in most browsers, but fails due to time-out error on certain networks and browsers.
The solution I think is to add a bounding box filter to this function. Could anybody tell me where to find the information or how to do this?

    var rootUrl = 'http://map.kecoviewer.com/geoserver/ows';

var defaultParameters = {
    service: 'WFS',
    version: '1.1.0',
    request: 'GetFeature',
    typeName: 'cresh:datazones_popup_mini',
    maxFeatures: 6505,
    outputFormat: 'text/javascript',
    format_options: 'callback: getJson',
    srsName:'EPSG:4326'

};

var parameters = L.Util.extend(defaultParameters);
$.ajax({
    jsonp : false,
    url: rootUrl + L.Util.getParamString(parameters),
    dataType: 'jsonp',
   jsonpCallback: 'getJson',
    success: handleJson
});

var group = new L.featureGroup().addTo(map);
var geojsonlayer;
var featureLayer = new L.GeoJSON();
var defaultStyle = {weight: 0,opacity: 0,fillOpacity: 0};
var highlightStyle = {color: '#f7f90e',weight: 3,opacity: 1,fillOpacity: 0};
var currentlayer;
function handleJson(data) {
    geojsonlayer=L.geoJson(data, {
        onEachFeature: function(feature, layer) {
            layer.setStyle(defaultStyle);

            (function(layer, properties) {
              layer.on("click", function (e) {
                if (currentlayer) currentlayer.setStyle(defaultStyle);
                layer.setStyle(highlightStyle);
                window.open("http:\/\/map.kecoviewer.com\/birtviewer\/run?__report=report\/datazone_story_ig.rptdesign&Datazone="+properties.dz,"story")
                currentlayer = layer;
            });

            })(layer, feature.properties);
}
    }).addTo(group);

}

function getJson(data) {
console.log("callback function fired");
}


texts: Normalizing raster using raster calculator raises TypeError

texts:
In ArcGIS 10.1 I'm trying to normalize a section of a DEM (digital elevation model) from original values to values ranging from 0 to 100. First I clipped the raster DEM with an irregulary shaped polygon.
Now I'm tying to change all the values to a range of 0 to 100 for further analysis.

I'm using the Raster calculator tool with this equation:

("DEM" - min("DEM")) * 100 / (max("DEM") - min("DEM")) + 0

I keep getting TypeError: 'Raster' object not iterable. This is the log I get:

   ERROR 000539: Error running expression: rcexec()
    Traceback (most recent call last):
      File "<expression>", line 1, in <module>
      File "<string>", line 5, in rcexec
    TypeError: 'Raster' object is not iterable

I also tried running the command on the unclipped raster, unsuccessfully.


texts:
I've got 2 rasters , its seems that 1 has the data equalized (all values from 1-65535 are there) and the 2nd don't .

Is there a process I can run in order to equalize the data in the 2nd raster ?


texts:
I will be offline (for mapping purposes) for a few weeks, but would like to continue working on my projects. Is there a way I can develop a leafletjs app while offline? I don’t want to create an offline map storage as I really don’t want to look at the maps. I just want leafletjs to not throw any errors while I code other (non-map) parts of the app.

update: to be clear, I want to be able to work with markers, coordinates, zoom, pan, etc. Just don't have to actually view the map tiles as there would be no connection to retrieve the maps.


texts: Stuck while Using GeoServer WPS with QGIS

texts: I have been using Geoserver WFS and WMS with Qgis for a longtime. Now i am trying to integrate WPS... I downloaded and activated WPS plugin for Geoserver. And i succesfully got result of buffer example with WPS Request Builder of Geoserver. But when i try with QGIS WPS-Client i never got result. It sends the data to server but never comes back the result. it stucks... I tried official 52North WPS server and i got result succesfully. So what should i do ???


texts: External geojson not showing in OpenLayers 2?

texts:
I'm new to OpenLayers. I tried some recipes found here to load an external geojson file as an overlay on a map I'm creating and it's not working. The map as well as my base layer are loading, but no overlay...

Can anyone help?

Both files, trains and gares.geojson are in the same folders as my HTML.

This is an extract from gares.geojson (in case the problem's there):

{
"type": "FeatureCollection",

"features": [
{ "type": "Feature",
"properties": { "NO_STATION": 101.0, "NOM": "Gare Lucien-L'Allier", "NoAbord08": 18, "Parking": 0 },
"geometry": { "type": "Point", "coordinates": [ -73.570806, 45.494761 ] } }
,
{ "type": "Feature",
"properties": { "NO_STATION": 214.0, "NOM": "Gare Delson", "NoAbord08": 41, "Parking": 1 },
"geometry": { "type": "Point", "coordinates": [ -73.540803, 45.366832 ] } }

]
}

I literally copy paste those two recipes and only changed the data source.


texts: How do I prevent labeling outside of the current atlas feature polygon?

texts: I want to create an atlas which is based on region polygons. I use inverted polygons and rule based styles (inverted polygons) for the polygon layer with the rule  $id  = @atlas_featureid which works fine. I have also a point layer and the features of the point layer are not shown due to the inverted polygons but their labels are shown. I want to prevent the labels outside of the current atlas feature polygon. Which expression should I use for the rule-based-labeling? I found this post (Filtering out labels surrounding active feature in Atlas) and tried the suggested rule-based filtering expressions but it didn't work for me.


texts: Adding points or hover tips from CSV file

texts:
I am using
QGIS version
3.4.1-Madeira
and want to add either a point on a map or hovertips on data that is in the attribute table.

What is the best method to do this?


texts: Why is QGIS Field calculator button greyed out?

texts:
I'm using QGIS to attempt to calculate the area enclosed by a set of points that I have capture used a GPS device.

I have imported the coordinates from a CSV but when trying to open the field calculator from the layer's attribute table the button is greyed out and cannot be clicked. Here is a screenshot of what I am seeing

Any ideas?


texts:
I would like to know how to transform the following code in something that only uses PyQGIS. I haven't found an example showing how to create raster files with PyQGIS yet.

I'm making a processing algorithm in QGIS 3. The variable input is a QgsRasterLayer.

#Getting input attributes
band_count = input.bandCount()
rows = input.height()
cols = input.width()
pixelType = input.dataProvider().dataType(1)

#Getting the destination file name and file format
destfile = self.parameterAsOutputLayer(parameters, 'OUTPUT', context)
fileformat = GdalUtils.getFormatShortNameFromFilename(destfile)

The actual raster creation is made below. This is the part I want to change to use only PyQGIS instead of GDAL. The variable classified is a numpy array made with scikit-learn.

driver = gdal.GetDriverByName(fileformat)
#Number of bands should be 1
band_count = 1
outRaster = driver.Create(destfile, cols, rows, band_count, pixelType)
outRaster.SetGeoTransform(dataset.GetGeoTransform())
outRaster.SetProjection(dataset.GetProjectionRef())
outRaster.GetRasterBand(1).WriteArray(classified)


texts: Using Query Builder with Joined table in QGIS?

texts:
When I join a layer with another layer all the attributes are correctly displayed in the first layer.

However when I use the Query builder (Layer - Query) and then use Zoom to layer, nothing is displayed.

I have added a vector layer with attribute values from Excel and joined it with another layer consisting of polygons. The polygon layer then displays all the attributes in its table from the Excel imported layer. When I  want to view only selected polygons using the Layer > Query using the joined attribute values and use Zoom to Extent nothing happens.

Rule-based styling seems to be the answer.


texts: How to add predefined values in exact order for value map field in PyQGIS and add empty field?

texts:
I have a problem with predefined values and their order in 'value map' field. Wnen I tried this How to create predefined values in a field in PyQGIS? it works, but I'm not able to keep order (even when I use OrderedDict instead of dict).
The second problem is how to insert also empty value.

    val_map = OrderedDict()
    for value,label in sorted(val['values'].items()):
        val_map[label.decode("utf-8")] = label.decode("utf-8")
    print(val_map)
    vlayer.setEditorWidgetV2( index, 'ValueMap' )
    vlayer.setEditorWidgetV2Config(index, val_map)

result:   OrderedDict([(u'silnice I.t\u0159', u'silnice I.t\u0159'), (u'silnice II.t\u0159', u'silnice II.t\u0159'), (u'silnice III.t\u0159', u'silnice III.t\u0159'), (u'MK I.t\u0159', u'MK I.t\u0159'), (u'MK II.t\u0159', u'MK II.t\u0159'), (u'MK III.t\u0159', u'MK III.t\u0159'), (u'MK IV.t\u0159', u'MK IV.t\u0159'), (u'UK', u'UK')])

But you can see, that there is alphabetical order in value map


texts: How to preserve topology while editing polygons in OpenLayers?

texts:
I’m facing this issue: I have a chloropleth map based on a set of polygons with some data. Using WFS-T (postgis, geoserver, openlayers) I am able to edit the geometry of the polygons. This works fine.

My polygons look like state borders.  They have shared border vertices. When I edit the geometry of some polygon (e.g., move vertices using modifyFeature), I need to affect the geometry of any neighboring polygon. There shouldn’t be holes or overlaps.

Is it possible to do it?

Thanks