texts: Automating creation of ServiceDefinition from Python script tool of ArcPy without broken project data source?

texts:
Summary:
I can create a ServiceDefinition through ArcCatalog or ArcMap, but when I try to generate the ServiceDefinition through python, it does not copy the /data folder with it, thus the script fails when published to the server.

Details:

I'm using ArcGIS 10.2.2.3552

My toolbox script has one string input parameter and one string output parameter

I have a standard toolbox with custom Python script tools that I need to put on an ArcGIS server. The requirements from the system administrator is that I give him a .sd (ServiceDefinition) file.

I can successfully create the .sd file through ArcCatalog by running my python tool and then sharing it as a service. The problem with this is that it is a very long process, and every step (run the tool, right click result share as service, analyze the tool, stage the tool) takes a very long time.

What I am trying do is create a python script that will fully automate the process of creating a ServiceDefinition file.

When I AnalyzeForSD, I keep getting an error in the analyzeMessages saying:

  ERROR 001238: Script myTool contains broken project data source: data

I can't figure out how to resolve this error.

I figured I could just ignore the error and try to create the ServiceDefinition anyway. So then the .sd file is created, but it doesn't copy /data folder.

I have tried copy_data_to_server=True and =False, but neither one actually copies the data.

Any ideas?

Here's my directory structure

Toolbox.tbx
generateServiceDefinition.py
ToolData

myTool.py
data

MapData.mxd
TextData.txt

Here is the script for creating the ServiceDefinition:

import arcpy
import sys
import os
scriptPath = sys.path[0]
toolboxPath = os.path.join(scriptPath, 'Toolbox.tbx')
sdDraftPath = os.path.join(scriptPath, 'Toolbox.sddraft')
sdPath = os.path.join(scriptPath, 'Toolbox.sd')

#import the toolbox
arcpy.ImportToolbox(toolboxPath, 'customModule')
#execute the tool
result = arcpy.myTool_customModule('hello')

try:
    os.remove(sdDraftPath)
except:
    pass
try:
    os.remove(sdPath)
except:
    pass

arcpy.CreateGPSDDraft(result, sdDraftPath, 'Toolbox',
                        server_type="ARCGIS_SERVER",
                        copy_data_to_server=False,
                        folder_name=None,
                        summary='Summary', tags='Some Tags')

analyzeMessages = arcpy.mapping.AnalyzeForSD(sdDraftPath)

if analyzeMessages['errors'] != {}:
    # if the sddraft analysis contained errors, display them
    print analyzeMessages['errors']
    raise Exception('There were errors while generating the service definition')

arcpy.StageService_server(sdDraftPath, sdPath)

Here is myTool.py. Obviously this isn't my actual script, but it is an example that reproduces the problem.

import arcpy
import os
import sys

input = arcpy.GetParameterAsText(0)
#print input to the console
arcpy.AddMessage('You passed in ' + input)
scriptPath = sys.path[0]
#make a variable pointing to my data so that ESRi's code will realize
#that this data needs to be packaged with the script
dataPath = os.path.join(scriptPath, 'data')
arcpy.AddMessage(dataPath)
mxdPath = os.path.join(dataPath, 'MapData.mxd')
#read a file from the data folder
dataFilePath = os.path.join(dataPath, 'TextData.txt')

with open(dataFilePath, 'r') as fin:
    arcpy.AddMessage(fin.read())

mxd = arcpy.mapping.MapDocument(mxdPath)
arcpy.AddMessage(arcpy.mapping.ListLayers(mxd))
#write something to the output parameter
arcpy.SetParameterAsText(1, arcpy.mapping.ListLayers(mxd))


texts: How to duplicate polygons in QGIS 3.0 with python

texts:
My background is in arcpy and Esri so all of the Qgis is new to me.

So I cannot figure out how to duplicate polygon with python. I want to draw one polygon and then be able to use it to make more and make a grid of them like I am showing below in the picture.

What should I do?
What resources should I look at? (I have looked at the documentation but doesn't make sense to me)

from qgis.core import *

layer = QgsVectorLayer(r"path/Test_file.shp", "Test_file", "ogr")

feats = [ feat for feat in layer.getFeatures() ]

The_Features = layer.getFeatures()
num=0
for i in The_Features:
    geom = i.geometry()

    Point_1 = geom.vertexAt(0)
    raw1 = Point_1.asWkt()
    Point_1_x = float(raw1[7:29])
    Point_1_y = float(raw1[31:51])
    print(raw1)

    Point_2 = geom.vertexAt(1)
    raw2 = Point_2.asWkt()
    Point_2_x = float(raw2[7:29])
    Point_2_y = float(raw2[31:51])
    print(raw2)

    Point_3 = geom.vertexAt(2)
    raw3 = Point_3.asWkt()
    Point_3_x = float(raw3[7:29])
    Point_3_y = float(raw3[31:51])
    print(raw3)

    Point_4 = geom.vertexAt(3)
    raw4 = Point_4.asWkt()
    Point_4_x = float(raw4[7:29])
    Point_4_y = float(raw4[31:51])
    print(raw4)
    print(Point_4_x)
    print(Point_4_y)
    num+=1
    break

 print(num)
 # Set the provider to accept the data source
 prov = layer.dataProvider()
 points = [QgsPointXY(Point_1_x,Point_1_y), QgsPointXY(Point_2_x,Point_2_y),
 QgsPointXY(Point_3_x,Point_3_y), QgsPointXY(Point_4_x,Point_4_y)]

 # Add a new feature and assign the geometry
 feat = QgsFeature()
 feat.setGeometry(QgsGeometry.fromPolygonXY( [points]))
 prov.addFeatures([feat])

 # Update extent of the layer
 layer.updateExtents()


texts: Problem with shooting star with pgRouting: doesn't respect the one ways

texts:
I used osm2po to create a table for pgRouting.
I notice a problem: for a route with dijkstra or a* the route is correct and respect the one ways. But with shooting*, it doesn't respect the one ways.

The route is from 11 baulacre geneve (46.2159082, 6.1394522) to coq inde cologny (46.2068106, 6.1919172).

Here are the queries:
dijkstra/a*:

SELECT topo.id, topo.osm_name, topo.geom_way
FROM Shortest_Path('SELECT * FROM ch_2po_4pgr', 263951, 271357, true, true)
AS path, ch_2po_4pgr AS topo
WHERE path.edge_id = topo.id;

shooting*:

SELECT topo.id, topo.osm_name, topo.geom_way
FROM Shortest_Path_shooting_star(
   'SELECT *, geom_way as rule, reverse_cost as to_cost FROM ch_2po_4pgr', 321089, 309056, true, true)
AS path, ch_2po_4pgr AS topo
WHERE path.edge_id = topo.id;

Did you notice the same problem with shooting* in others situations?

Here is a picture.
The orange line is for a*.
The red line is for shooting*.
Since with shooting* with use edge_id instead vertex_id, maybe the algo starts at the beginning of the way (starting point is in pink), but the problem with the one ways is in the pink circle: the red line uses the "rue des Alpes" which is a one way.

Thanks to all for your help!


texts: Guide to aerial or space data types and availability

texts:
I am overwhelmed by the variety of types of satellite and aerial data (images?) available online. For example, at the USGS Earth Explorer data page, I see AVHRR; CEOS;Ikonos and OrbView; EO-1; GLS; HCMM; LPDACC (which contains AirMOSS, ASTER, MODIS, SRTM, ...)

At different times, I am seeking data for different purposes. Sometimes I want clear, high-resolution images; other times it may be something else.

Is there a good guide to the aerial and satellite data products available for download? I would like to read about what data types are good for what mapping purposes, and then see references to tutorials or instructions on how to process/use the specific data type.

I use QGIS.


texts: How to draw a polygon with preset length in OpenLayers?

texts: For a project I am working on, I need to be able to put a polygon on a map using OpenLayers, based on a variable point that can be set by a user. The polygon is an isosceles triangle of which two legs should be exactly 1 kilometer and should originate from a point set by the user under an angle of 90 degrees.  Does anyone know how this can be done?


texts: Change/Hide the pen width for line shape file based on rule in QGIS

texts:
I have line shapefile, in which I have a condition to display line's width based on a rule as shown in the image below. But in this case I was not able to define line width to 0, it is taking hairline as minimum by default. Is there any other way to hide the other features in print-composer, so that I can print only selected features?

Basically I need to print only the green line with hidden black lines:


texts: Saving Split Layer segments to separate Shape files in QGIS?

texts:
I might be missing something obvious here.  Using the Split Features tool I've split a line feature into 2 parts, one of which is now showing as being selected (yellow highlight).

How do I save each part to a new Shapefile?


texts: Georeferencing with Quantum GIS and reprojecting to EPSG:900913

texts:
I'm trying to reproject a map to EPSG:900913
The image I'm reprojecting is the following http://i.imgur.com/HmfLW.png

I'm using the cities as reference points obtaining the coordinates from Google Maps.

I'm entering the coordinates for each of the reference point this way:

So X is longitude in degrees and Y is latitude in degrees.

I'm using the following transformation settings and target srs:

Given these settings the "Generate GDAL script" action gives me the following script:

gdal_translate -of GTiff -gcp 39.084 83.0534 12.9639 55.6032 -gcp 569.16 509.313 24.0161 49.838 -gcp 109.924 489.771 14.3701 50.1065 -gcp 420.153 334.656 20.9839 52.2413 -gcp 426.26 581.374 21.2256 48.7345 -gcp 549.618 142.901 23.8623 54.9019 miasta.bmp miasta_translated.bmp
gdalwarp -t_srs epsg:900913 -r near -order 1 -co COMPRESS=NONE  miasta_translated.bmp miasta_modified.tif

It is strange it is not adding the -t_srs parameter to the gdalwarp command, so I've added it myself.

I expect the miasta_modified.tif to contain the Google Maps projection of the original map but instead it contains the following which doesn't look like the Google Maps projection the most noticeable difference being that it is compressed in the vertical direction: http://i.imgur.com/9DMyT.png (I've converted it to PNG for the image hosting to accept the file)

What am I doing wrong?


texts: Click on an object to open an image using a QGIS action

texts:
I'm working on this answer from @Nathan W to find how to deal with missing files, as they return a error message if QGIS doesn't find any.

The initial workflow is to be able to click on an object (here a point which is a borehole) and open the related crosscut, that I have renamed from the borehole name. If the crosscut doesn't exist, I'd like to display a message saying this crosscut is not referenced. All the images are in jpg format and are in a folder named "Fichiers_Coupes_SCAN" ad my image name is in the attribute table in a column named "ID_SONDAGE".

My code is the following :

from os import startfile
from os.path import abspath, dirname, join, isfile
from PyQt4.QtGui import QMessageBox
proj = QgsProject.instance()
urfile = str(proj.fileName())
path = join(abspath(dirname(urfile)),"Fichiers_Coupes_SCAN", "[% "ID_SONDAGE" %]"".jpg")
import os.path
if os.path.isfile(path) :
    startfile(path)
else:
    QMessageBox.information(None, "Error : ","Can't find the crosscut")

I'd like the error message to display the name of the borehole, but I must do something wrong and I get error messages when inserting the variables.


texts: Shapefile Polygon does not display in QGIS

texts:
I am attempting to import a Shapefile into QGIS.  When I import it, the boundary is shown on the list, but when I try to view it in the mapping area, it does not show up.

I can get it to show the nodes if I go into edit mode, so I know the data is there.  I've have the same problem in v2.8.2 and v2.12.3.  I can run some of the polygon tools from the tool box such as SAGA Polygon properties, and it will then show the polygon when the new layer is created.

I am able to load these same polygons into another program and they work properly.  I've saved the files to my Google Drive and it can be linked below.

Anyone have any idea what I might be doing wrong.  I've gotten other shapefiles to show in QGIS with no real issues previously.  I'm very new to QGIS.

link to Boundary data Shape Data


texts:
I've developed a simple web map using GeoServer and OpenLayers in WMS format. The styles are created in Geoserver using SLD and send image tiles to the client. In order to display a point layer, I've used PointSymbolizer and external graphics which shown properly on the map. The problem is the there is not any VendorOptions to control icon overlapping and some points close to each other are drawn on top of each other as shown in the picture.

The following is the of SLD for the point icons. How can make icons to not overlap with each other (the SpaceAround is not the option).

<se:Rule>
  <se:Name>fuel</se:Name>
  <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
   <ogc:PropertyIsEqualTo>
        <ogc:PropertyName>nearby_category_fk</ogc:PropertyName>
        <ogc:Literal>103</ogc:Literal>
      </ogc:PropertyIsEqualTo>
  </ogc:Filter>
  <se:MaxScaleDenominator>136494</se:MaxScaleDenominator>
  <se:PointSymbolizer>
    <se:Graphic>
        <!--Plain SVG fallback, no parameters-->
        <se:ExternalGraphic>
            <se:OnlineResource xlink:href="poi_new_cats/fuel-day.png" xlink:type="simple"/>
            <se:Format>image/png</se:Format>
        </se:ExternalGraphic>
        <!--Well known marker fallback-->
        <se:Mark>
            <se:WellKnownName>square</se:WellKnownName>
            <se:Fill>
                <se:SvgParameter name="fill">#000000</se:SvgParameter>
            </se:Fill>
            <se:Stroke>
                <se:SvgParameter name="stroke">#ffffff</se:SvgParameter>
                <se:SvgParameter name="stroke-width">0.5</se:SvgParameter>
            </se:Stroke>
        </se:Mark>
        <se:Size>20</se:Size>
    </se:Graphic>
   <se:VendorOption name="labelObstacle">true</se:VendorOption>
  </se:PointSymbolizer>
</se:Rule>


texts: Directly Rendering raster to Google Map API from Google Earth Engine

texts:
I am using Google Earth Engine Python API for calculating NDVI and then sending the NDVI result to Google cloud storage. After that I am displaying the result (NDVI tiles stored in Google bucket/cloud storage) at front end using Google Map API. But exporting the result from Earth Engine to bucket takes average 5 to 10 minuets depending on my ROI. Now Earth Engine API itself can render the output within few seconds.

My question, is it possible to directly render the raster output on Google Map API from Google Earth Engine (without storing the result to bucket.. as storing result in Google bucket takes much more time)?


texts: How accessing a specific feature and comparing it with the next one?

texts:
Let's imagine to have a vector layer whose attribute table is made of an ID field, a field "A" (string) and a field "B" (float).

I am trying to access a specific feature (e.g. [0]) in field A, comparing the selected value with the one after([1]) in the same column and, based on this, doing something.

If I were working on a spreadsheet, I'd write: IF(A0=A1; B0=x;"").
Based on the final hint enclosed in a previous discussion (link), which is quite closed, I tried to use it but:

it moves among fields (horizontally) and not features (vertically)
and I should activate it for each feature (but what about if there are many features?)

SCRIPT TRIAL:

features = layer.getFeatures()
attrs = feat.attributes()
attr = attrs[0]
Index=2 # field index
request = QgsFeatureRequest().setFilterFid(0)
#feature index
feat = layer.getFeatures(request).next()
print feat.attributes()[Index]
#thus, the result returned is the first feature in third field
Index=Index+1
#I change the index to the 1st.one on the right
print feat.attributes()[Index]
#thus, I obtain the 1st feature of the 4th field


texts: NVIZ not working with QGIS 2.16 and GRASS 7 (mac)

texts:
I am trying to display 3d data in NVIZ (ideally i would like both raster and vector files in 3d) but when I run the NVIZ tool, nothing is returned. Im running QGIS 2.16.3 with the GRASS 7 plugin installed (GRASS mapset is activated).

I am using projected rasters stored in QGIS.

**I notice that when i put these rasters into a GRASS location they are not recognized by the NVIZ tool at all.

I have read that some people have had success installing a stand alone version of GRASS, but I would rather avoid this if possible.


texts: How to determine if a .obj 3D Mesh file is geo-referenced or not?

texts: I have a .obj file output from OpenDroneMap. But I cannot determine if it is geo-referenced or not. How can I determine that?


texts:
I have used rasterize function from GDAL in order to rasterize a shapefile that I have according to existing raster.
The problem is that I'm not able to find or to show the resulting raster.
This is the rasterize function I am using:
def rasterise_me(raster, vector, attribute,
                fname_out="", format="MEM"):
    """Rasterises a vector dataset by attribute to match a given
    raster dataset. This functions allows for the raster and vector
    to have different projections, and will ensure that the output
    is consistent with the input raster.

    By default, it returns a handle to an open GDAL dataset that you
    can e.g. `ReadAsArray`. If you want to generate a  GTiff on disk,
    set format to `GTiff` and `fname_out` to a sensible filename.

    Parameters
    ----------
    raster: str
        The raster filaname used as input. It will not be overwritten.
    vector: str
        The vector filename
    attribute: str
        The attribute that you want to rasterize. Ideally, this is
        numeric.
    fname_out: str, optional
        The output filename.
    format: str, optional
        The output file format, such as GTiff, or whatever else GDAL
        understands
    """
    # Open input raster file. Need to do this to figure out
    # extent, projection & resolution.
    g = gdal.Open(raster)
    geoT = g.GetGeoTransform()
    nx, ny = g.RasterXSize, g.RasterYSize
    srs = g.GetProjection()
    min_x = min(geoT[0], geoT[0]+nx*geoT[1])
    max_x = max(geoT[0], geoT[0]+nx*geoT[1])
    min_y = min(geoT[3], geoT[3] + geoT[-1]*ny)
    max_y = max(geoT[3], geoT[3] + geoT[-1]*ny)
    # Reproject vector to match raster file
    vector_tmp = gdal.VectorTranslate("", vector, format="Memory",
                                    dstSRS=srs)
    # Do the magic
    ds_dst= gdal.Rasterize(fname_out, vector_tmp, attribute=attribute,
                        outputSRS=srs, xRes=geoT[1], yRes=geoT[-1],
                        outputBounds=[min_x, min_y, max_x, max_y],
                        format=format, outputType=gdal.GDT_Int32)
    return ds_dst

then I apply it on my data:
rasterize=rasterise_me(r"the/path/to/my/raster.tif", r"the/path/to/my/shape.shp", "class",fname_out='Rasterize.tif')

I couldn't find the results so I had tried to define the path also for fname_out like this:
rasterize=rasterise_me(r"the/path/to/my/raster.tif", r"the/path/to/my/shape.shp", "class",fname_out='this/is/the/path/for/my/future/raster/Rasterize.tif')

but still couldn't find it.
My endgoal is to be able to find my lost rasterize result raster.


texts: How to cut the end of line in a QGIS style?

texts:
I want to create a style in QGIS, in which one line is offset from the other, but their ends form one point (as in the picture). How do can you remove the ends of the green line?

In QGIS 3.2 there is an EXTEND function that extends the ends, but does not shorten them...


texts:
QGIS 1.8 will not start on Windows 7

Basically a similar issue as the above post. Version 1.8 fails to install with the standalone installer; with OSGeo4w, while it doesn't run into trouble during the installation, ultimately fails as well (same error message as below when opening the app).

Standalone installer has the following problems:

At the end of the installation, it reaches "qgis-common.bat" execution and stalls; meanwhile, Task Manager shows "find string" GREP utility running in the background and repeatedly crashing every minute or so. Ending the process results it immediately starting back up.

If the installation is forced to quit or completed, trying to open Qgis results in the error message that qwt5.dll is missing.

QGREP crash signature:

Problem Event Name:APPCRASH
Application Name:   findstr.exe
Application Version:    6.1.7601.17514
Application Timestamp:  4ce78c0a
Fault Module Name:  findstr.exe
Fault Module Version:   6.1.7601.17514
Fault Module Timestamp: 4ce78c0a
Exception Code: c0000005
Exception Offset:   00003c76
OS Version: 6.1.7601.2.1.0.768.3
Locale ID:  1033
Additional Information 1:   0a9e
Additional Information 2:   0a9e372d3b4ad19135b953a78882e789
Additional Information 3:   0a9e
Additional Information 4:   0a9e372d3b4ad19135b953a78882e789


texts: Creating shapefile showing footprints of Rasters?

texts:
I have around 1,000 satellite images in tiff format, and I want to create a shapefile which will serve as an index to the rasters. This is something similar to a raster catalog, but I do not want to build a raster catalog.

Some obstacles I can forsee, is that the image are georefrenced, so they are not rectangular in shape (I am talking about the data area).

To clarify, I require the polygon to cover only the non-zero (or non nodata) pixels of the raster, and not the entire rectangular raster. Most of the answers so far, give a rectangular polygon, which covers the data, as well as the non-data pixels.

My Image

Result given by  tools I have examined (like raster catlog, various Arcscripts, custom Python script given in one of the answers):

Result that I want:


texts: In QGIS my expression aggregation in attribute forms slow works

texts:
There are two layers: Parcels and Parking.
In parking has attribute: parcels_objects, this is lists of parcels intersect of parking.
When I create a parking object, I want to calculate all parcels containing parking lots and write a list of their numbers to the parcels_objects parking field attribute.
To do this, in the attribute forms for the parcels_objects field, I created an aggregate expression in default value:
aggregate(
    layer:='Parcels',
    aggregate:='concatenate',
    expression:="parcel_number",
    concatenator:=', ',
    filter:=intersects($geometry,geometry(@parent))
)

This expression works, but it is very slow because (6 seconds before opening the attribute form) there are 50,000 objects in the parcel layer.
Spatial index was created in both tables:
CREATE INDEX parcels_sx ON public.parcels USING gist (geom);
CREATE INDEX parking_sx ON public.parking USING gist (geom);

Dividing workflows into 2 different ones will speed up the assignment of the parcels_objects attribute, for example:

create objects first
call the process of assigning parcel numbers

But the idea is that the business process of entering information was continuous, for example:

the graphic object is created
all the attribute data for the current object are filled in (some attributes are filled in automatically) with the possibility of correction.

The virtual field in my case is not entirely suitable, since I want to calculate the parcels_objects with each edit.
Form opening slowdown occurs when the option is checked: Apply default value on update:

How can you make this expression run faster?


texts: geocoding addresses specifically into NAD1983 or converting X,Y into 1983 after geocoding

texts:
I'm going bananas trying to geocode addresses into a NAD1983.

I have ArcMap 10.1 and I'm using the default address locator (North American 10.0 from arcgis online).  No matter how low I set the sensitivity, it gives me "There was an error trying to process this table" when I specifically choose NAD1983.  If I allow it to create the X,Y coordinates in WGS1984, then I get coordinates...BUT my other coordinates (I'm going to measure distance/proximity) are in NAD1983.

Any suggestions? I can't figure out how to change them from WGS1984 to NAD1983 when they're in a layer, although I have changed the data frame successfully it seems...

Any guidance is very appreciated!


texts: Transferring Z attributes to polygon layer in QGIS 3.10

texts:
I have a point layer with the z attribute (and there is a field with the same z value). I have digitised a polygon, snapping to the point layer. How can I transfer the z value from the point layer to each polygon vertex?

Edit

They're a bunch of points surveyed with a total station of the outline of a building, so I have connected the dots to make the actual footprint polygon, but each polygon vertex has a 0 Z attribute, and I want to transfer the Z attribute from each point to each corresponding vertex.


texts: Get list of all tables in CartoDB account with SQL statement?

texts:
How can I get a list of all available tables in my own authenticated CartoDB account with a SQL statement?

This returns null:

SELECT table_name as t FROM information_schema.tables WHERE table_schema = 'public';

Alternatively, could it be done with CartoDB.js?


texts:
I recall having a label-toggle button in the toolbar in earlier QGIS (2.x). If the selected layer had labels configured, then clicking this button would toggle their visibility off or on. The labels' formatting, placement etc would remain unchanged.. they'd just go visible or invisible. It was very useful functionality, particularly when I had to zoom out of one area and move to another. I'd toggle the labels off when I didn't need them - it also sped up the tool's response time. Can't figure out how to get that button again in QGIS 3. I can't remember if it was some specific plugin that did it, and now can't find any such plugin.

What do I need to do to get the toggle labels button back?

EDIT: Same question asked and answered here: Pause Labeling in QGIS?

And yes it was the Deactivate/Active Labels plugin. Missing it! It's currently not 3.x-compatible. Started a request on their repo.


texts:
I have a DEM .tif data. I am trying to convert it with gdal2xyz to ASCII data.

The result is always wrong. X and Y values are correct but Z value is wrong.

On the next image below, you can see the parameters. I changed even the band but it did not help. Any suggestions?

I am working with EPSG: 31468 | DHDN / 3-degree Gauss-Kruger zone 4.

You can find the result (CSV-data) in the next image:


texts:
I have a model running on the processing framework it produces a serie of vector layers as output(alt-0, alt-1, alt-2, alt-3, and so on). I have prepared the model keeping the order of these alternatives from zero to 5.
But when running the model  these vector layer are listed pretty much randomly (see pic below). I tried to re order the .model with a notepad editor but nothing changes.
Is there a way to set the order on the window of the model for the outputs?


texts: Creating parcels of equal areas in QGIS?

texts:
I want split a polygon of approximately 100ha (land) into parcels of 500m2 and I'm   lost.

Is there a tool, extension, or code that you can share with me to do it using QGIS?


texts:
This should be simple: have a shapefile in UTM z16 DatumNad27 and want to transform it to Datum WGS84 latlon
I am using Qgis 2.6

Sorry I am going to clarify, with the steps taken:

1- XY data of a landfarm, I converted to shape and load it in Qgis. (NAD27_UTM 16N)

2- I have recorded gpx data, a pathway to a volcano (amazing by the way). Load gpx data in Qgis.(WGS84)

3- When loading google satellite or OSM Landscape (WGS 84/Pseudo Mercator), gpx data aligns very well. But line derived from XY is not aligned.

4- Of course I am using "transformation on the fly". But it does not aligned.

5- When projecting XY line, using SAGA algorithm,and specifying WGS84 as targetCRS, it produces no changes (see image)


texts:
I've got a PostGIS databases, with lots of schemas. I can easily create a new user that will allow them to update and delete features. However when I get the user to create a new feature I get the following error in QGIS:-

  PostGIS error while adding features: ERROR:  new row for relation "conservation_areas" violates check constraint "enforce_geotype_geom

I've used the following SQL to GRANT the user (anne) privilages that will allow her to update and delete, but I guess its the insert bit that is letting me down.

GRANT ALL ON SCHEMA conservation TO anne;
GRANT SELECT ON ALL TABLES IN SCHEMA conservation TO anne;
GRANT UPDATE ON conservation.conservation_areas_gid_seq To anne;
GRANT UPDATE ON conservation.conservation_areas TO anne;

The postgres user has no such issues when edited data, therefore I'm missing something somewhere when i'm granting the user permissions. Is there a simple SQL that will allow a user to do everything that the Postgres user can do but to a specific tables (and not all tables) within a schema?

Med


texts: ArcSDE edits not showing up in ArcGIS Server Geodata Service?

texts:
I've got an ArcGIS Server GeoData service against an SDE MSSql database.  When changes are made to the source data, the edits don't appear to show up in the GeoData service.

If I restart the service, edits are available via the geodata service.

Is there something else I should be doing to allow edits to always be available?

Further investigation revealed that the service did reflect changes, though how long you have to wait is variable.  Sometimes the delay is two minutes, sometimes 30 minutes.

At this stage, to force the GeoDataService not the cache any data, I've set the "Virtual Output Directory" of the service to be none.  No caching, no problems.  Though not ideal, as when no edits are happening, we don't get the performance improvements that caching offers.


texts:
I am looking to generate points along a series of connected lines at regular intervals, but I wish to keep this interval even when a line is branching into 2 or 3 other lines. All the regular tools I tested in ArcMap or QGIS generate points on a feature basis meaning that the interval restart at each node.

On the picture I drew road network (blue lines), direction flow (orange arrows) and desired points (red dots). The first point is at 50m from the beginning of the first line, then a second 50m further but then there is a node within the next 50m. The usual tools would then start "fresh" on the two branches but I want to keep the 50m distance between the last point before the node and the first point on each branch. How could this be done? I suppose some topology may be needed...

I am looking for a free solution, ideally with Python, QGIS, SAGA, GRASS, GDAL or alike. In last resort, I could use ArcMap.


texts: Interpolate shape point data with integer data into raster

texts:
I have a shapefile I brought into QGIS with a lot of point data.  One of the data layers represents an integer value that is assigned to each position.

I am attempting to create a raster map of the data but I don't want the interpolation process to calculate non-integer values.  I simply want it to fill in the spaces as "1", "2", etc. and change basically where they meet.

Eventually I want to create a Contour map that shows these different zones.  The only way I know how to do that currently is to create a raster first and then convert.  If there is another way to do that, I'd be happy to try that instead.

Is there any way to automate this process or do I have to draw these polygons manually?


texts: Using the CSS style format in GeoServer, how can you sort-by on a value returned from a function, i.e. the x value of a polygon’s centroid

texts:
Work in progress
I’m currently working on a CSS style in GeoServer where the US states are extruded based on population (PERSONS). In order to get a 2.5D effect a combination of the isometric and offset functions needs to work in conjunction and drawn in a specific order.

Dificulty
From the picture above one can see that NY (New York) are drawn on top of PA (Pennsylvania). I therefore would like to specify the drawing order based on the x value of the polygon’s centroid.   I’m also all ears if there is a better solution because I don’t believe sorting by the x value will eliminate all problem areas.
Tried and Failed
To get the x value of the polygon’s centroid, one can use this function:
[getX(centroid(the_geom))].
When applying this to the sort-by function I get the following error:
sort-by: [getX(centroid(the_geom))];

Cannot turn this value into a literal: Expression [expression=getX([centroid([the_geom])])]
Full code
@mode "Flat";

* {
  stroke: #000000;
  stroke-opacity: 1;
  stroke-width:3;
  z-index:0;
}

* {
  fill: #7B7B7B;
  fill-geometry: [isometric(the_geom, PERSONS/8000000)];

  stroke: #636363;
  stroke-geometry: [isometric(the_geom, PERSONS/8000000)];
  stroke-opacity:0.7;
  sort-by: PERSONS;
  z-index:0;
}

* {

  fill-geometry: [offset(the_geom, 0, PERSONS/8000000)];
  stroke-geometry: [offset(the_geom, 0, PERSONS/8000000)];
  fill: #CACCCD;
  stroke: #000000;
  stroke-opacity: 0.7;
  sort-by: PERSONS;
  z-index:0;
}

*{
  label-geometry: [(centroid(offset(the_geom, 0, PERSONS/8000000)))];
  label-anchor: 0.5 0.5;
  label: [STATE_ABBR];

  font-family: "Arial";
  font-fill:black;
  font-size:12;
  halo-color: white;
  halo-radius: 1.5;
}

So in this example I would like to replace the sort-by: Persons with sort-by: x value of polygon centroid.


texts: Basic unit testing for ArcGIS 10.0 Java Add-in

texts:
I've been assigned the task of translating an old ArcGIS VBA script into a Java add-in. I finished the rough translation of the code, now I'm starting the debug process. The first unit test I want to set up is the creation of an Application object. I'm getting the error "0x80040154 - Class not registered". Any ideas? I'm a very novice programmer so I'm open to any suggestions. Thank you.

Here's the code for my test:

@Test
public void testOnClick() {
    try {
        EngineInitializer.initializeEngine();
        AoInitialize ao = new AoInitialize();
        ao.initialize(esriLicenseProductCode.esriLicenseProductCodeArcInfo);
        IApplication app = new Application();
        assertNotNull(app);
    }
    catch (UnknownHostException e) {
        System.out.println("UnknownHostException: " + e.toString());
    }
    catch (IOException e) {
        System.out.println("IOException: " + e.toString());
    }
}


texts: Project point to map with Arcpy using GPS Coordinates

texts:
I am trying to create a new layer with just a single point using arcpy.  All that is provided is a gps lat and long.  It seems like it should be a simple problem but I just can't get it to work.  I've determined the error is that I'm not using the correct spatial reference (the point I get right now is next to the equator/meridian) but I'm not sure how to fix that.

import arcpy

arcpy.env.workspace = 'c:\\Scrap'
arcpy.env.overwriteOutput = True

latitude = 41
longitude = -71

p = arcpy.Point(longitude, latitude)
pg = arcpy.PointGeometry(p)
arcpy.CopyFeatures_management(pg, 'test.shp')


texts: Discrepancy between constructed line length and length recorded in attribute table in ArcGIS for Desktop?

texts:
Say that in an edit session I create a line of an exact length entered by keyboard - say 10'

If I later measure that line using the measure tool, snapping to each end of the line, it reads 9.99998'

Then when I check the Attribute Table immediately after completing the line, I find that the line length is 9.99998'.

This is an actual example I just tried, but it it seems to be consistent with all such measurements. They're all off - always less than what I entered as far as I can tell - by a small amount.

Now, I realize this is a very small difference, but shouldn't the editing tool build a line that is exactly 10' if that is the length I enter, and then shouldn't the measure tool "measure" the line to be 10' if that is the length I entered when creating the line?

What is going on here? And although for what I am likely to be doing it doesn't matter, is there a way to prevent this.

To clarify: I create the line in an ArcView edit session by clicking where I want the line to start, then right clicking and selecting "exact length" to type in the length of the line I want to create. I type in 10', and the next node is dropped, in theory at least, exactly 10' from the last point. Then I click F2 to complete the drawing. When I check the attribute table it shows that that the segment is 9.99998' rather than the 10' I specify. I first noticed when I was using the measuring tool - snapping at the beginning and the end of the line. I'm not exactly sure how I could clarify the more.

I am using ArcView 10.1 SP1
Data frame and data are both NAD_1927_StatePlane_New_York_Long_Island_FIPS_3104.
Data are created and stored in a File Geodatabase. Not sure what the tolerances are at the moment.


texts: AUTO:42003 Orthographic projection outside valid area

texts:
When running the following code to reproject a epsg:4326 countries file to orthographic:

newCRS = CRS.decode("AUTO:42003,1,-100,45");
mapFrame.getMapContent().getViewport().setCoordinateReferenceSystem(newCRS);

The following error occurs:
WARNING: Possible use of "Orthographic" projection outside its valid area.
Latitude 112°43.8'S is out of range (±90°).

Does anyone know how to get the orthographic projection to work?


texts: Providing user specified location for extracting packages in ArcGIS Pro?

texts:
In ArcMap the folder location used to unpack packages dragged into a map can be set using the Customize | ArcMap Options > Sharing tab:

I just looked for the equivalent option in ArcGIS Pro on the Project Page | Options but it has no Sharing tab:

Is the only way to control this by using the output folder of the Extract Package tool instead?


texts: How to merge street layers in QGIS for exporting to PostGIS?

texts:
Is there any way of merging lots of street data in QGIS and getting one layer for exporting it to Postgresql?  I know one method using mmqgis "Merge Layer" tool.  It works in Ubuntu but not in Windows. When I want to open my MapInfo .TAB file in Windows it opens but this time mmgqis tools do not work.  In Ubuntu it gives "*.TAB is not a valid or recognized data source "...

Can I solve this with some Python code or any other way?

Thanks


texts: How to serve web maps with user-specific datasets from ArcGIS Server

texts:
I am a new ArcGIS Server user and am trying to understand how our web mapping application specification would fit into the ArcGIS Server map services architecture.

Our application should support multiple users, each of whom will be working with a limited number of their own maps. Each map is made up of a common base map across all users, but with a user-maintained layer from a geodatabase that is specific to the individual user.

So should I be thinking in terms of a single map service which is consumed by all the users, and somehow for each request I swap in the correct layer for that user. Or rather should there be a separate map service for each user (and for each of the user's maps) - in which case I guess I would need a way to dynamically create & publish the map services to the server at runtime.

In either case I am looking at writing a Server Object Extension to provide the som of my custom application functionality.


texts: Geocoding API: house number precision for Europe (Germany/Austria), other than google?

texts:
I am looking for an affordable geocoding option with the following criteria:

commercial use
<500$/month
allows to store LatLng results in my DB
house Number precision for Austria and Germany

The only options with sufficient coverage seem to be BING and Google.
I have already received a quote from google for $11000/year for their cheapest commercial license and it would still not allow me to store the results.

Most other options that I am aware of are either not suitable for my area (e.g. mapquest) or are based on OSM which also lacks many house numbers where I need them.

It seems that I am pretty much out of luck and now i understand why google charges such ridiculous licensing fees... because they can.


texts: Spatial Index usage for points-in-polygon query in spatialite

texts:
I have two datasets, one of which contains millions of postcode points and the other which includes polygon data - county borders in this instance. I'd like to write an ogr2ogr command to export only the points which are contained by a specific polygon.

My source data is within separate SQLite databases, and my output needs to be a ShapeFile.

I've written a VRT file which produces the correct output:

<OGRVRTDataSource>
    <OGRVRTLayer name="pcodes">
        <SrcDataSource>MyPostcodes.sqlite</SrcDataSource>
        <SrcLayer>PostcodeData</SrcLayer>
    </OGRVRTLayer>
    <OGRVRTLayer name="polys">
        <SrcDataSource>Counties.sqlite</SrcDataSource>
        <SrcSql>
            SELECT GEOMETRY
            FROM CountyTable
            WHERE countyname="Nottinghamshire"
        </SrcSql>
    </OGRVRTLayer>
    <OGRVRTLayer name="res">
        <SrcDataSource>test.vrt</SrcDataSource>
        <SrcSql dialect="SQLITE">
            SELECT pcodes.GEOMETRY FROM pcodes, polys
            WHERE Within(pcodes.GEOMETRY, polys.GEOMETRY) = 1
        </SrcSql>
    </OGRVRTLayer>
</OGRVRTDataSource>

My ogr2ogr command is:

ogr2ogr -f "ESRI Shapefile" pointsinpoly.shp test.vrt -overwrite -sql "SELECT * FROM res"

Although this produces the correct ShapeFile, it takes a few minutes to run. I realise that I need to use the spatial indexes, but I'm stumped as to what the correct syntax is.

I have a SpatialIndex table on each of the referenced databases and I've confirmed they're both valid with select CheckSpatialIndex().


texts: Accounting for Colour Blindness when designing maps?

texts:
Depending on what you read, 7 to 15% of the general population has some form of colour blindness.

Can this be accounted for in  map designs and, if so, how?


texts:
I have multiple .tif files downloaded from Ortho Downloader for an Aerial View of Ward 21 in Saint Louis City, MO.  I need the .tif files to underlay my layer depicting the boundaries of Ward 21 as well as the parcels within Ward 21.

The .tif layers come projected to NAD_1983_2011_StatePlane_Missouri_East_FIPS_2401_Feet.  My parcel layer and Ward layer are both projected to NAD_1983_StatePlane_Missouri_East_FIPS_2401_Feet.  And yet when both are added to the same blank ArcMap fullscreen view shows them as projected at different locations - the .tif files on the right and the rest of my St Louis layers on the left, separated by several inches of space.  I have tried Geographic Transformations, Projecting each layer to 1983 again, and restarting my computer.  It will not let me project the .tif files (they simply don't show up in the Input options for Project).  I cannot figure out how to make these things lay on top of each other.  Please help!


texts: Tab widgets in QGIS plugin - run different scripts

texts: Is it possible to run separate python processes through tab widgets of the same QGIS plugin. I have created a QGIS plugin using the plugin builder. From QT designer I added two tab widgets to the plugin. I want to run different processes in both the tabs, one at a time though. Right now I can run the process in the first Tab, but I do not know how to add the script to the second tab and run it when the 'OK' button in second tab is clicked.


texts: KML file not showing anything on the google map

texts:
When i import this KML file: https://yadi.sk/d/G3DJywNAliflew

Nothing shows up on the map, just centering around Kazakhstan. This KML file has layer (like 2G cellular coverage map).

What i am doing wrong? Can i use this KML in a Yandex maps as well?


texts: Minimizing error when buffering unprojected data in ArcGIS Desktop?

texts:
If I buffer hundreds of unprojected WGS84 points across the US by say, 3 miles, will I get accurate areas returned (3 x pi x r) for each of the circles created?  I suspect not. If not, what would be some best practice to minimize error for each buffer.

Basically, I'd like to minimize error across the US (perhaps including Hawaii) and processing time (in ArcGIS, for instance).

For context, my goal is to extract demographics (also in unprojected geography) for each point lying within a constant distance of each point for the entire nation. My fear is that using unprojected data for this task means that if I'm comparing the demographic info captured by each buffer, I'm not comparing apples to apples.


texts: Leaflet marker cluster not working on GeoJSON with MultiPoint geometry

texts:
I have a GeoJSON file of fly tipping incidents that I have managed to display on a basic Leaflet map.

I am now trying to incorporate the marker cluster plugin, but I am having problems as it does not seem to be clustering the markers when zooming out.

My code is below, is there an obvious reason why it is not working? I am receiving no error messages.

(I followed the tutorial from http://maptimeboston.github.io/leaflet-intro/)

// initialize the map
  var map = L.map('map').setView([51.394845, -0.538516], 12);

  // load a tile layer
  L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png',
    {
      maxZoom: 17,
      minZoom: 9
    }).addTo(map);

  // load GeoJSON from an external file
  $.getJSON("data/incidents.geojson",function(data){
    // add GeoJSON layer to the map once the file is loaded
    var incidents = L.geoJson(data,{
      pointToLayer: function(feature,latlng){
        var marker = L.marker(latlng);
        marker.bindPopup('Fly tip reported on ' + feature.properties.date_rec + '
' + 'Reported waste type: ' + feature.properties.waste_type + '
' + 'Located on land type: ' + feature.properties.land_type);
        return marker;
      }
    });
    var clusters = L.markerClusterGroup();
    clusters.addLayer(incidents);
    map.addLayer(clusters);
  });


texts: I have preparing road symbology using Arcmap 10.1 version. When I was prepared attribute based   symbology as a layer shown on symbols were overlapped. I have attached an example image to show the problem.


texts: Getting started with mapping flow of natural gas using GIS?

texts:
I'd like to get involved with GIS, im a marine engineer and clueless about it and although i've tried to research it, I'm overwhelmed by the amount of programs available & information. I wanted to do a personal project, what i have in mind is try to map the flow of Natural Gas.

So how should i go about it to complete my project?

I need a free program to start; which one would you suggest?

Does it have forums / documentation support?

I'm alone on this, so the easier I can get info/answers the better. I'd rather get a local program ie. not web based. I run windows7.

As I understand it, I also need:

A world map with countries?
Some import/export data which I think I've found.
What else would I need to do the analysis from the program pov?
Some analysis tool to get the flows? Is that included in the program
or must I find an add-on or something like this?

My future targets include finding/plotting sea routes & also pipelines. Would the suggested program be capable of doing that? Would it automatically be able to plot the sea routes?

Sorry if I wasn't clear but what I'm asking is what free GIS program can do all of the above?

Is there a step by step guide to learn?

As suggested from the posts QGIS is a starting point & read its manual with a goal to do networks analysis. ArcGIS is the industry standard and when I have enough money I should consider it. Getting a book as well would help. Long-term doing python programming or something.


texts: Writing Python script for plugin to filter out high slope areas in QGIS

texts:
I am currently developing a Python plugin using the Python QGIS API that tries to automate basic constraint analysis. So far, I am able to add buffers to parcel, wetland, and floodzone areas. But, I would now also like to write some code that will identify high slope areas on a DEM layer. Doing this will allow me to further restrain my build area.

So far, I am able to save the DEM layer to a raster object like so -

rlayer = QgsProject.instance().mapLayersByName('38_77_3688')[0]

I have taken a look at the Python QGIS API documentation and have come across the interesting class class qgis.analysis.QgsSlopeFilter at https://qgis.org/pyqgis/master/analysis/QgsSlopeFilter.html. But I am having some trouble understanding how to use this class.

Can someone explain what QgsSlopeFilter does or point me to a more useful class that will allow me to filter out the high slope areas in the DEM layer?


texts: Plotting shapefiles with Pyshp in Python omits some shapes

texts:
I want to plot a cloropleth map of Scotland with Pyshp, but only a fraction of the shapes are plotted. The shapefile of Scotland has 32 shapes, and it is available here. The values dictating the colors are in the following dictionary:

scotland=
{'Aberdeen City': 16.52,
 'Aberdeenshire': 13.805,
 'Angus': 12.252,
 'Argyll and Bute': 11.158,
 'Clackmannanshire': 8.346,
 'Comhairle nan Eilean Siar (Western Isles)': 9.999,
 'Dumfries and Galloway': 7.72,
 'Dundee City': 8.377,
 'East Ayrshire': 9.844,
 'East Dunbartonshire': 7.741,
 'East Lothian': 10.661,
 'East Renfrewshire': 10.334,
 'Edinburgh (City of)': 11.593,
 'Falkirk': 11.534,
 'Fife': 10.318,
 'Glasgow City': 10.65,
 'Highland': 11.066,
 'Inverclyde': 11.145,
 'Midlothian': 10.58,
 'Moray': 10.697,
 'North Ayrshire': 9.173,
 'North Lanarkshire': 10.755,
 'Orkney Islands': 10.618,
 'Perth and Kinross': 10.425,
 'Renfrewshire': 10.138,
 'Scottish Borders': 8.231,
 'Shetland Islands': 10.181,
 'South Ayrshire': 7.863,
 'South Lanarkshire': 10.846,
 'Stirling': 10.021,
 'West Dunbartonshire': 8.272,
 'West Lothian': 10.488}

The script I use is faulty:

import shapefile
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Polygon
from descartes.patch import PolygonPatch

sf_sc = shapefile.Reader("C:\\MyName\\MyFolder\\Scotland\\Scotland_with_council_areas.shp")

#Scotland
plt.figure()
ax = plt.axes()
ax.set_aspect('equal')

#Fieldname of interest: Hygiene
fld_sc = sf_sc.fields[5:] #That's where Hygiene is stored in the shapefile attribute table
field_names_sc = [field[0] for field in fld_sc]
fld_ndx=4

#Setting global min/max of the records
maxrec=-9999
minrec=1e21
for rec in sf_sc.iterRecords():
    maxrec=np.max((maxrec,rec[fld_ndx]))
    minrec=np.min((minrec,rec[fld_ndx]))
maxrec=maxrec/2.0 #Fixing the upper saturation limit

for shapeRec in sf_sc.iterShapeRecords():

    #Getting shape geometry and records
    shape=shapeRec.shape
    rec = shapeRec.record

    #Selecting polygon colors (RGB values) based on the Hygiene values
    R = 1
    G = (rec[fld_ndx]-minrec)/(maxrec-minrec)
    G = G * (G<=1) + 1.0 * (G>1.0)
    B = 0

    print [R,G,B]

    #Number of parts in each polygon (we never know how the shapefiles were digitised)
    nparts = len(shape.parts)
    if nparts == 1: #Shape with a unique part
       polygon = Polygon(shape.points)
       patch = PolygonPatch(polygon, facecolor=[R,G,B], edgecolor=[0,0,0], alpha=1.0, zorder=2)
       ax.add_patch(patch)
    else: #Looping over the parts of each shape and plot separately
       for ip in range(nparts):
           i0=shape.parts[ip]
           if ip < nparts-1:
              i1 = shape.parts[ip+1]-1
           else:
              i1 = len(shape.points)

              # build the polygon and add it to plot
              polygon = Polygon(shape.points[i0:i1+1])
              patch = PolygonPatch(polygon, facecolor=[R,G,B], alpha=1.0, zorder=2)
              ax.add_patch(patch)

plt.xlim(5000,500000)
plt.ylim(500000,1300000)
plt.show()

It is faulty because for some reason, only 6 of the 32 shapes are plotted as a result:

While we all  know that the shape of Scotland is quite different:

Printing the [R,G,B] sequence returns:

[1, 0.11049362536776752, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 0.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 0.016345210853220233, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 0.42759071592023518, 0]
[1, 0.39490029421379469, 0]
[1, 1.0, 0]
[1, 0.48447204968944063, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 0.50866296175220671, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]
[1, 1.0, 0]

What am I missing?


texts:
First off, these questions are relevant:

https://dba.stackexchange.com/questions/7900/pros-and-cons-of-having-one-table-per-entity-vs-having-a-shared-table-with-an-e/7907#7907

https://dba.stackexchange.com/questions/55951/deciding-when-to-extend-to-child-tables

I'm quite new to database administration so I'm looking for some advice on how to set up a geographic workorder system. Specifically, I have several tables, each with a geographic feature type.

A workorder will consist of a description, date, etc, a feature type, and may be related to one or many features of one type. Each feature may be in many workorders, so I realize this is a many-to-many relationship, which confuses me even more.

Ultimately this data will be used in a python based ORM (SQLAlchemy). Flask will handle the creating and modification of the workorders. A mapping server will be accessing the data also, to map the workorders to their related features geometry.

The features are created in an Esri geodatabase and dumped into PostgreSQL regularly (using ogr2ogr).

Here is what I've come up with so far:

Table Feature x - Each table of features with feature info and geometry
Table workorder - holds the workorder info, date, description...
Table workorder_feature - holds the workorder_id and feature_id

Workorder will have one to many relationship with Workorder_Feature, and since the Feature tables are getting droppped/recreated regularly I'll just have to recreate views that join the feature to their relevant workorder.

Another option I thought about was to have a table of workorders for each feature. But this doesn't seem like the best way to solve this though, if anyone has any alternatives I'd love to hear them.


texts:
I've created a plugin icon which can toggle and is connected to the run() function. At the moment, whenever a user clicks on the icon, it executes the run() function.

How could I make it so that it only executes this function when its toggle state isChecked()?

A simple example:

def add_action(self,icon_path,text,callback,checkable=False,enabled_flag=True,add_to_menu=True,add_to_toolbar=True,status_tip=None,whats_this=None,parent=None):
    icon = QIcon(icon_path)
    action = QAction(icon, text, parent)
    action.triggered.connect(callback)
    action.setEnabled(enabled_flag)
    action.setCheckable(checkable)

    if status_tip is not None:
        action.setStatusTip(status_tip)
    if whats_this is not None:
        action.setWhatsThis(whats_this)
    if add_to_toolbar:
        self.toolbar.addAction(action)
    if add_to_menu:
        self.iface.addPluginToMenu(
            self.menu,
            action)

    self.actions.append(action)
    return action

def initGui(self):
    icon_path = ':/plugins/Example/icon.png'
    self.add_action(
        icon_path,
        text=self.tr(u'Example'),
        callback=self.run,
        checkable=True,
        parent=self.iface.mainWindow())

def run(self):
    print 'Plugin loaded.'
    #if checkable == True:
    #   print 'Checked'
    #else:
    #   print 'Unchecked'


texts: How to connect to CouchBase with ogr2ogr?

texts:
I am trying to load a Shapefile to CouchBase (version: 2.5.0 DEV edition, installed on Ubuntu), but with no success.
According to this documentation it is possible like this:

  ogr2ogr -f couchdb
  "couchdb:http://some_account.some_couchdb_server.com" shapefile.shp

so I tried like this:

  ogr2ogr -f couchdb "couchdb:http://localhost:8092" shapefile.shp

or

  ogr2ogr -f couchdb "couchdb:http://localhost:8092/osm" shapefile.shp

(where osm is an existing databucket)
and also with permissions:

  ogr2ogr -lco UPDATE_PERMISSIONS=ALL -f couchdb
  couchdb:http://localhost:8092 shapefile.shp

then I tried simply to connect:

  ogrinfo -ro "couchdb:http://localhost:8092"

which gives me the error:

  FAILURE:  Unable to open datasource 'couchdb:http://localhost:8092'
  with the following drivers. [list of drivers including couchdb]

I also tried 127.0.0.1 instead of localhost and also with port 8091 (and every combination of these).

Can anyone verify that the couchdb ogr2ogr driver is working? And if yes, what am I doing wrong?
Maybe the problem is that it is CouchBase instead of CouchDB?
 My GDAL version is 2.0.0dev
I am grateful for any suggestions.


texts: How to get images displayed on ArcGIS Online Popups

texts: I have added images to show up in an html popup on ArcMap, but when I bring it in to ArcGIS Online, they don't appear. Is there a way to get the images to show up there? Currently on ArcGIS Online the option offered only shows the same image for all the points.


texts: False DescribeFeature Request URL in GeoExt QueryPanel

texts:
I got the GeoExt QueryPanel working so far, it recognizes the layers available on my Geoserver but it is not able to get the attribute columns for the dataset. I configured an osgeo proxy on my Apache Server, I use the override-ext-ajax.js and have the OpenLayers.ProxyHost = "/cgi-bin/proxy.cgi?url="; class configured.

A DescribeFeatureType request is sent when selecting the layer I want to filter. The same happens when trying to choose the attribute column in the query panel, but there are none displayed.

The sent request URLs look like this:

http://localhost/cgi-bin/proxy.cgi?url=
http%3A%2F%2Flocalhost%2Fcgibin%2Fproxy.cgi%2Fwfs%3F
VERSION%3D1.1.0%26REQUEST%3DDescribeFeatureType%26TYPENAME%3Dcases%26query%3D

The strange part here is, that the request is sent via the proxy.cgi twice. The request URL should look like this, going passed geoserver´s WFS:

http://localhost/cgi-bin/proxy.cgi?url=
http%3A%2F%2Flocalhost%3A8080%2Fgeoserver%2Fwfs%3F
VERSION%3D1.1.0%26REQUEST%3DDescribeFeatureType%26TYPENAME%3Dcases%26query%3D

Does anybody have a clue what´s going wrong here and why the request URL is wrong? The strange thing is that even if I remove the Open.Layers.ProxyHost class and the override-ext-ajax.js from my application, the request is sent that way. This leads to the conclusion, that the source of error can be found in the QueryPanel.js file, yet I couldn´t spot it.

Can it, on the other hand be related to my server?

Other requests sent in the application (GET & POST) work fine. The (relevant) code looks like this, maybe you can spot the error.

OpenLayers.ProxyHost = "/cgi-bin/proxy.cgi?url=";

var owsUrl = "/geoserver/wfs";

var symbolizer = OpenLayers.Feature.Vector.style["default"];

symbolizer.graphicName = "square";

var store;

Ext.onReady(function() {

            store = new GeoExt.data.WFSCapabilitiesStore({
            url: Ext.urlAppend(owsUrl,
               "VERSION=1.0.0&REQUEST=GetCapabilities"),
            protocolOptions: {version: "1.1.0"},
            autoLoad: true,
            listeners: {
                load: initApp
            }
            });
        });

function initApp() {
        var vecLayer = new OpenLayers.Layer.Vector("Query", {
            eventListeners: {
                featuresadded: function() {
                    var extent = vecLayer.getDataExtent()
                    mapPanel.map.zoomToExtent(extent)
                }
            }
        });

 var queryPanel = new gxp.QueryPanel({
            title: "Query",
            region: "center",
            width: 350,
            bodyStyle: "padding: 10px",
            autoScroll: true,
            layerStore: store,
            bbar: ["->", {
                text: "Query",
                handler: function() {
                    queryPanel.query();
                }
            }],
            listeners: {
                query: function(panel, store) {
                    var rule = new OpenLayers.Rule({
                        filter: panel.getFilter()
                    });
                    rule.symbolizer[queryPanel.symbolizerType] = symbolizer;
                    var style = new OpenLayers.Style("", {rules: [rule]});
                    wms.mergeNewParams({
                        sld_body: new OpenLayers.Format.SLD().write({
                            namedLayers: [{
                                name: panel.selectedLayer.get("name"),
                                userStyles: [style]
                            }]
                        })
                    });

                    store.bind(vecLayer, {
                        initDir: GeoExt.data.FeatureStore.STORE_TO_LAYER
                    });
                }
            }
        });


texts:
What's the difference between using CRS defined as

+proj=longlat +datum=WGS84 +no_defs

and a CRS defined as

+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

Will using one or the other cause any problems. In the first case, how can GIS project properly without knowing which ellipsoid the data are referred to?


texts: Creating ALOS PALSAR 1.0 level JAXA Orthorectified K & C Mosaic interferometry

texts:
JAXA is providing free access to PALSAR 50m Orthorectified K & C Mosaic.

The files are not in the typical PALSAR format. I'm looking for apps that can read and export this data.

Is it possible to create interferometry using these products?

They formatted that product differently than raw imagery, with a _HDR file and an _IMG file for each band.
http://www.eorc.jaxa.jp/ALOS/en/kc_mosaic/kc_map_50.htm


texts: Print Username in print composer of QGIS

texts: Is it possible to use any expression to get the username of machine/computer and path of project while printing the documents in print composer of QGIS.


texts: GeoDjango - save two Point fields, one for geography and one for geometry

texts:
Most of the explanations I've read for geography and geometry seem to suggest that one must only use either geometry or geography, and that the only way to use both is to cast/convert between one and the other on-the-fly. But is this actually the case?

Why not save two points?

class Point(models.Model):
    googlemaps_placeid = models.CharField(max_length=254)
    geography = models.PointField(geography=True)
    geometry= models.PointField(srid=3348) #Statistics Canada Lambert projection

Most of the queries from my users will be of the type "Show me all points within roughly 20km from this point, and order them from nearest to furthest". Even though the data has points all over Canada, the user is only interested in points within a small locality with any given query, so the geometry PointField would suffice.

If distance measurements in meters are required (as opposed to simply ordering them from nearest to furthest), then the geography PointField will be used.

Are there any flaws with this plan? If not, what is the actual syntax for saving the two PointFields? (I've searched all over for this since GeoDjango docs are severely lacking, and all I've found is this with no accepted answer: https://stackoverflow.com/questions/24921545/how-to-assign-to-a-django-pointfield-model-attribute)


texts: need Matching Algorithm which takes graph Topology into account

texts:
Problem:

I have a problem (illustrated below) with matching two datasets.

On the left I have a GIS containing drain network data (in blue), including diameter, length and topology information (from- and to manholes).

I have received some video drain survey data which has neither coordinates nor a reliable key to the original GIS data (illustrated on the right in green). It does however contain useful attribute information like diameter and length and topology information. Not all drains have been surveyed.

I need to match the two datasets. I can do this by matching similar diameters and length attributes (providing a score) and pairing up drains using bipartite best matching. I do not however use the topology which could be extremely valuable information and ensure the matching occurs in the correct order. My approach could therefore deliver poor results (imagine a network where lots of drains have similar lengths and diameters).

I would like to know if anyone knows of an algorithm (Python implementation?) applicable to this type of problem.

I am considering generating the score by also comparing up- and downstream pipes but this just seems like an extension of a bad idea...


texts: QGIS2THREEJS Installation Windows Problem

texts: I just install qgis2threejs on the plugin (i use qgis 3.4.3 madeira) but I can't activate and I get caution about this. I try to reinstall qgis but it didn't work. after i click ok then qgis close automatically. anyone can help me?


texts: How to disassociate roughness from slope

texts:
How to measure topographic roughness irrespective of slope in ArcGIS? I have used many methods to generate roughness layers, including TRI, VRM (which is supposed to be less correlated with slope), slope aspect index, slope variance, slope standard deviation etc.

So far, all of the output layers are extremely similar to each other and to my slope layer. They do not seem to be clearly differentiating between areas of steep and rugged terrain and those which are steep but smooth.

How can I generate a roughness layer which is not so dependent on slope?


texts:
I m currently trying to build app in the Intramaps roam. IntraMaps Roam is a simple data collection application built using QGIS. Roam is built as a standalone fully packaged Python application with everything included. QGIS is not required to be installed on the users machines. see this url for more information.

At the moment I m trying to add a topo background layer. This backgroundlayer (brtachtergrondkaart) is a WMTS provided by the dutch gouverment and is open source. The url to this WMTS layer
is: https://geodata.nationaalgeoregister.nl/tiles/service/wmts?request=GetCapabilities&service=WMTS

I have added the WMTS layer to the QGIS project that the app uses. When i start the app and open the map I get the following message: Failed to load brtachtergrondkaart.

I can't find in the documentation of Intramaps roam if a WMTS layer is supported. Does anybody know if this is integrated or do I have to make this feature in the program?


texts: Ubuntu 14.04 and QGIS 2.8.1 unmet dependency of libgdal .so.1-1.11.1

texts:
I know that QGIS and unmet dependencies seem to be a common theme on this site, but this most recent one is a known issue. Through the standard Ubuntu updates in 14.04, at some point QGIS 2.8.1 has an unmet dependency with libgdal .so.1-1.11.1 and gets autoremoved. The solution on the ubuntu-gis list is to upgrade to the nightly 2.9.x version.

However this is not without its problems: Strange circles have appeared on QGIS map? The circles appear for me in 2.9.x and I can't get them to go away.

There is currently no accepted answer to that question, so I wonder if there is some other way to install this vital package.


texts: Drawing lines from point to polygon at angle to shortest line between same point and polygon?

texts:
I have a polygon with a point outside. I was able to calculate the shortest distance between the point and polygon and draw a line for the same. Now, I need to draw radiating lines from the point, on either side of the shortest line, such that they meet the polygon. The shortest line is the reference and the radial lines are decided by angle between them. The angles are  theta,2*theta,3*theta etc

I didn't come up with a code for whatever I have managed so far; rather, I used a bunch of functions using arcpy. I'm still in the initial process of understanding the possibilities and what part of my requirement needs custom coding.

The steps that I followed are as follows:

Used Near_analysis to find the nearest point on the polygon from my point. It added the location(X and Y) of the nearest point in the attribute table
Using the coordinates of the nearest point and the coordinates of my point, I implemented the XYtoLine_management to draw the line.

At this point I'm stuck. I've been trying bearing to distance but that uses geographic north as reference.

It not much of a code as it has several errors which are irrelevant, with a few results generated as per my needs.

I'm still exploring some options.


texts: GeoExt and Mapfish print on other domain: can't get how to proxy

texts:
I have OpenLayers ProxyHost configured

OpenLayers.ProxyHost= "form/proxy?url=";

, the url config of Mapfish Print Provider is configured as http://path/to/mapfish/print, and the doc tells that I can do that as it is accessible throug the proxy

This property requires that the print service is at the same origin as the application (or accessible via proxy).

, but the lib still tries to go to other origin directly and fails.

I cannot find where to configure GeoExt or ExtJS to go there through proxy.

The Mapfish Print doc tells to set some parameter of a mapfish variable, I did that but obviously with no success.

Please, help!

(Actually, it's GeoExt 2 for ExtJS 4, where the PrintProvider seems being  renamed as GeoExt.data.MapfishPrintProvider)


texts: Does it matter whether I use ArcCatalog or Access to compact a geodatabase?

texts:
I've been tasked with compacting some of the personal geodatabases (MDB files) in our office every few weeks. These reside on a couple of separate file servers.

I have been using ArcCatalog to do this task, but it is very slow to launch, and switching from directory to directory also causes the title bar to display "Not Responding" for a while. All in all, the process is painfully slooooow...

However, I remember performing some compact tests on these files using Microsoft Access several months ago, and it didn't have these same problems. And the databases seemed to behave normally afterwards. However, I decided to adopt ArcCatalog because I assumed that it might be better suited to compacting these special kind of MDB files.

Was this a correct assumption? Does it matter whether I use ArcCatalog vs. Microsoft Access when I want to compact my geodatabases? Is the compact process any different between the two?


texts: Adding features in the layer with reference to another layer

texts:
I have built a tool for creating line in the exact center of two selected lines through this code:

import math

def mainCl():
    ftr_geom = []
    select = iface.activeLayer().selectedFeatures()
    for ftr in select:
        geom = ftr.geometry()
        if geom.type() == QGis.Line:
            f_geom = geom.asPolyline()
            ftr_geom.append(f_geom)

l1x1=ftr_geom[0][0][0]
l1y1=ftr_geom[0][0][1]
l1x2=ftr_geom[0][1][0]
l1y2=ftr_geom[0][1][1]
l2x1=ftr_geom[1][0][0]
l2y1=ftr_geom[1][0][1]
l2x2=ftr_geom[1][1][0]
l2y2=ftr_geom[1][1][1]
pr = iface.activeLayer().dataProvider()

d1=math.sqrt(math.pow((l2x1-l1x1),2)+math.pow((l2y1-l1y1),2))
d2=math.sqrt(math.pow((l2x2-l1x1),2)+math.pow((l2y2-l1y1),2))
d3=math.sqrt(math.pow((l2x1-l1x2),2)+math.pow((l2y1-l1y2),2))
d4=math.sqrt(math.pow((l2x2-l1x2),2)+math.pow((l2y2-l1y2),2))

if d1<d2:
    case1()

else:
    case2()

def mid(pt1, pt2):
   x = (pt1.x() + pt2.x())/2
   y = (pt1.y() + pt2.y())/2
   return QgsPoint(x,y)

def case1():
    ftr_geom = []
    select = iface.activeLayer().selectedFeatures()
    for ftr in select:
        geom = ftr.geometry()
        if geom.type() == QGis.Line:
            f_geom = geom.asPolyline()
            ftr_geom.append(f_geom)

l1x1=ftr_geom[0][0][0]
l1y1=ftr_geom[0][0][1]
l1x2=ftr_geom[0][1][0]
l1y2=ftr_geom[0][1][1]
l2x1=ftr_geom[1][0][0]
l2y1=ftr_geom[1][0][1]
l2x2=ftr_geom[1][1][0]
l2y2=ftr_geom[1][1][1]
pr = iface.activeLayer().dataProvider()

midpt1= mid(QgsPoint(l1x1,l1y1), QgsPoint(l2x1,l2y1))
midpt2= mid(QgsPoint(l1x2,l1y2), QgsPoint(l2x2,l2y2))
nwl = QgsFeature()
nwl.setGeometry(QgsGeometry.fromPolyline([midpt1,
midpt2]))
pr.addFeatures( [ nwl ] )

def case2():
    ftr_geom = []
    select = iface.activeLayer().selectedFeatures()
    for ftr in select:
        geom = ftr.geometry()
        if geom.type() == QGis.Line:
            f_geom = geom.asPolyline()
            ftr_geom.append(f_geom)

l1x1=ftr_geom[0][0][0]
l1y1=ftr_geom[0][0][1]
l1x2=ftr_geom[0][1][0]
l1y2=ftr_geom[0][1][1]
l2x1=ftr_geom[1][0][0]
l2y1=ftr_geom[1][0][1]
l2x2=ftr_geom[1][1][0]
l2y2=ftr_geom[1][1][1]
pr = iface.activeLayer().dataProvider()

midpt1= mid(QgsPoint(l1x1,l1y1), QgsPoint(l2x2,l2y2))
midpt2= mid(QgsPoint(l1x2,l1y2), QgsPoint(l2x1,l2y1))
nwl = QgsFeature()
nwl.setGeometry(QgsGeometry.fromPolyline([midpt1,
midpt2]))
pr.addFeatures( [ nwl ] )

mainCl()

Now I want to modify my tool..I want that I select a layer and select its two lines..The new centerline should be drawn in the layer of my choice not in the same layer..To choose the layer in which I want to draw line is through this select box:

And the layers it loads should be the map layers I added in the layer box..
Overall task will be that I select two lines in one active layer and the centerline will be drawn in the the layer I selected in that box.

Can anybody help me how to perform this task?


texts: Obtain attributes of a point layer that falls within a polygon layer of different CRS

texts:
I have created my own vector grid polygons in QGIS 2.12.2. This creates a new layer to my project, but it's projected in a "generated projection".
I want to view the attributes of points that fall within the generated grids.

I have tried to save it as an ESRI shapefile under WGS84 however it disappears from my map once it is saved.

Is there a way to re-project those grids into WGS84?


texts: How to create a suitable filtered .osm for osm2po using OpenStreetMap Overpass API?

texts:
I successfully created a (.osm) set suitable for routing in a desktop GIS over a large custom area thus keeping only primary network (way object) using the Overpass API (see below).

I'm trying this approach because my study area is crossing national and local borders, so other approach to parse a full planet file might take a lot of resources... which I don't have...

However the resulting OSM XML was rejected by osm2po because the file tr_raw.2po is missing after the first reading phase (see log below).

Is there a way to produce through Overpass API a suitable file (but not too heavy!) to use along with osm2po?

Or can we modify osm2po behavior to work without relations so no error is thrown?

Thanks and sorry for my poor English ;)

Here is the overpass request:

<osm-script>
  <query type="way">
    <has-kv k="highway" regv="junction|motorway|motorway_junction|motorway_link|primary|primary_link|secondary|secondary_link|tertiary|tertiary_link|trunk|trunk_link"/>
    <bbox-query e="8.1024169921875" n="48.34529727896014" s="47.331377157798244" w="6.8060302734375"/>
  </query>
  <union>
    <item/>
    <recurse type="down"/>
  </union>
  <print/>
</osm-script>

And here is the log:

LOG    Opened at Thu Apr 25 11:59:14 CEST 2013
INFO   Running osm2po 4.7.7 with cmd=tjsp
       SingleTileMode Uncompressed - 1 359M
INFO   Starting Tiler at Thu Apr 25 11:59:14 CEST 2013
INFO   Reading from
       mypath\extract.osm
INFO   Using parser de.cm.osm2po.converter.OsmXmlParser
INFO   166 865 of 166 865 nodes extracted - 1 354M
INFO   0 of 0 relations extracted - 1 354M
DEBUG  Writing TileFile:tm_info.2po
DEBUG  Prefix:extract Ways:0 Relations:0
       Tiles[1] SingleTileMode
       [0] Id:0 (s090-w180) Coords:166 520 Complex:345
INFO   Tiler finished at Thu Apr 25 11:59:16 CEST 2013
INFO   Starting Joiner at Thu Apr 25 11:59:16 CEST 2013
DEBUG  Reading TileFile tm_info.2po
DEBUG  Prefix:extract Ways:0 Relations:0
       Tiles[1] SingleTileMode
       [0] Id:0 (s090-w180) Coords:166 520 Complex:345
DEBUG  Capacity for Shared Ids is 3
INFO   Caching relations from tr_raw.2po - 1 358M
FATAL  Exception at Thu Apr 25 11:59:16 CEST 2013
       Message: java.io.FileNotFoundException: extract\tr_raw.2po (Le fichier spécifié est introuvable)
       Class: java.lang.RuntimeException
       Stacktrace:
       de.cm.osm2po.primitives.InStreamDisk.<init>(InStreamDisk.java:30)
       de.cm.osm2po.converter.Joiner.cacheOsmRelations(Joiner.java:106)
       de.cm.osm2po.converter.Joiner.run(Joiner.java:171)
       de.cm.osm2po.Main.runJoiner(Main.java:191)
       de.cm.osm2po.Main.main(Main.java:126)
INFO   ------------------------------------------------------------------

       For more help enter
       java -jar osm2po-core-4.7.7-signed.jar -h

LOG    Closed at Thu Apr 25 11:59:16 CEST 2013


texts: Averaging elevation for fishnet points in ArcGIS for Desktop?

texts:
I have a problem with adding average elevations for fishnet points. First of all I surveyed an area boundary of 16km2 using total stations. I got nearly 800 points to finish the boundary in this way. I have created fishnet points with the total station extent points as the limit for fishnet points.

I want to add averaging elevations for the the fishnet points with that total station points as base points. Both are in shapefile, so am not looking for suggestions involving raster processes.

There are 3000 fishnet points. Without any idea I am adding elevations for those 3000 points by averaging using calculator.


texts: Can QGIS geopackage files with multiple layers exist as a single file with the .gpkg extension?

texts:
My question is about the basic structure of a .gpkg file. When I make a .gpkg file, they only ever contain one layer, and I cannot add additional layers to them. However, sometimes when I have added and then exported the vector layers as .qgz files and the completed map as .qgt files, all left in the same directory as the .qpkg file, they will show up in the layers or associations of the .gpkg file, IF I have saved the .gpkg file after performing all these operations.
As an example, following the above to create these files and saving them all in the same directory, would show two layers when I opened the .gpkg file:
X.gpkg - contains polygons [layer 1]
Y.gqz - contains points [layer 2]
Y.gqt - contains themes for map made [not seen as layer, but accessible via SHOW LAYOUT MANAGER]
Is there a way that I could save all this so that it was just in one .gpkg file? Or instead is the .gpkg file always found with its corresponding layer files as separate files?
I've watched numerous videos on YouTube, read the documentation but am still fundamentally confused about this.


texts: Remove username and password from saved arcgis layer (.lyr) files?

texts:
How can one remove the saved username and password from ArcGIS layer (.lyr) files?

A scripted solution preferred but even a point and click method will be welcome.


texts: QGIS create link line features between table points by point names


texts:
I have a multiband raster (3 bands) that I want to display in QGIS.  I can display each band as a singleband without trouble, but if I try display from multi-band (properties in style) I get very little of the image.  Tips for what I am doing wrong?  I am using the option to load the max/min values from extent.  See screen shot of the properties selected below, as well as a screen shot of just the blue band

ETA: Checking that this really is multi-band and not single band.  I have copied the properties using the check that @ahmadhanb has shown in their response.  3 bands are mentioned.

Driver
GDAL provider
GTiff
GeoTIFF
Dataset Description
filepath/SRTM_W84.tif
AREA_OR_POINT=Area
Band 1
RepresentationType=ATHEMATIC
STATISTICS_COVARIANCES=11535.08822294172,8768.945936995617,6876.945116165363
STATISTICS_MAXIMUM=255
STATISTICS_MEAN=195.17337427807
STATISTICS_MINIMUM=0
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV=107.40152787237
Band 2
RepresentationType=ATHEMATIC
STATISTICS_COVARIANCES=8768.945936995617,8122.704384270356,6278.396582608941
STATISTICS_MAXIMUM=255
STATISTICS_MEAN=209.64313615503
STATISTICS_MINIMUM=0
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV
0
Data Type
Byte - Eight bit unsigned integer
Pyramid overviews
Layer Spatial Reference System
+proj=longlat +datum=WGS84 +no_defs
Layer Extent (layer original source projection)
115.0000000000000000,4.9990704361995046 : 130.0001196658648723,20.0000000000000000
Band
Band 1
Band No
1
No Stats
No stats collected yet
Band
Band 2
Band No
2
No Stats
No stats collected yet
Band
Band 3
Band No
3
No Stats
No stats collected yet =90.126047202073
Band 3
RepresentationType=ATHEMATIC
STATISTICS_COVARIANCES=6876.945116165363,6278.396582608941,5652.621835039329
STATISTICS_MAXIMUM=255
STATISTICS_MEAN=169.71246008137
STATISTICS_MINIMUM=0
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV=75.183920056348
Dimensions
X: 18521 Y: 18522 Bands: 3
X : 9261,Y 9261
X : 4631,Y 4631
X : 2316,Y 2316
X : 1158,Y 1158
X : 579,Y 579
X : 290,Y 290
X : 145,Y 145
Origin
115,20
Pixel Size
0.000809898,-0.000809898
No Data Value 0 Data Type Byte - Eight bit unsigned integer
Pyramid  overviews Layer Spatial Reference System
+proj=longlat +datum=WGS84 +no_defs
Layer Extent (layer original source projection)
15.0000000000000000,4.9990704361995046 :
0.0001196658648723,20.0000000000000000
Band Band 1 Band No 1 No Stats No stats collected yet
Band Band 2 Band No 2 No Stats No stats collected yet
Band Band 3 Band No 3 No Stats No stats collected yet


texts:
I import a dted file into a postgis enabled database with the -R option.
I also gave the complete directory path to the file.
I also have the raster and gdal related environment variables on the system.
When executing a sql query that must use the dted data I get the following error:

rt_band_load_offline_data Cannot open offline raster

With the full path of where the file is.
Any ideas what the problem could be, or how to resolve this?


texts: How can I subtract one (very very large) set of features from another?

texts:
I'm trying to subtract a very large set of features (water areas) from another, relatively small set of features (land area).

In principle, I should be able to do this in QGIS using Vector > Geoprocessing Tools > Difference, but as the very large set of water features is very large (a 1.7GB shapefile), this is taking... a while. (Given that the progress bar on running a small subset of it has not budged in 14 hours, I am reasonably convinced it has crashed.)

What's the best way of doing this? Can I do it with gdal or ogr or a python library? Will it be more likely to work if I merge all the separate water features into a single feature? (Currently there are about 2.2 million, and I don't actually care about the water features, only the land.)


texts: I am looking for a way to do fast KNN queries on a dataset of 10+ million 3D points (and growing). I was hoping PostGIS would do the trick with its new <-> and <#> operators, but they don't seem to work in three dimensions. Does anyone know any workarounds that would extend the KNN GiST index to 3D? Right now I'm using a "normal" spatial index with the &&& operator, but that creates all the issues that led to the addition of the KNN index in the first place.


texts: Points within certain distance to each other

texts:
I've got a layer with more than 400.000 Points.
Now I want to select those within a distance of 100m to each other, to
verify that points in such a distance don't have different IDs.

How can I do this?


texts:
I have seen Where to get school zone shapefiles? which seems to be the closest to my need.

When I download the docs, either I have no ways to filter the document or I don't see all the schools. I have been trying for a week now.

I can download individual files here, but it definitely does not have all schools.

Where can I obtain all school boundaries (starting with CA)?

As I wrote in the very first paragraph, there is a question very similar to this one, but it does not answer my question for two reasons:

the Census offers School District and also some Elementary School boundaries, but they are not accurate down to each individual school.
I downloaded a huge file from that question, but it has no description (some generic name for each file), it is for the whole country, and I don't have any software that allows me to see the file. I need to click on every school in the country until I find my local school, and chances are these two sources are actually just one, and it will not have the individual elementary, middle and high school.

Image 1 is the city school district boundaries (which I also need, so that is helpftul).

Image 2 has some Elementary schools, not for all cities, and most importantly, it does not match reality. Each city has many elementary schools, and the map brings just one for "some" cities.


texts: Setting a string Output parameter from the python script tool in ArcGIS

texts:
import sys, arcpy
arcpy.env.workspace = "C:\Scratch"
Environment_Index = arcpy.GetParameterAsText(0)
Measure = arcpy.GetParameterAsText(1)
if Environment_Index == "Percentage of Environmental Hazard":
    if Measure == "0-25%":
        Expression = '"per_Envi_H" >0 AND "per_Envi_H" <25'
    elif Measure == "25-50%":
        Expression = '"per_Envi_H" >25 AND "per_Envi_H" <50'
    elif Measure == "50-75%":
        Expression = '"per_Envi_H" >50 AND "per_Envi_H" <75'
    elif Measure == "75-100%":
        Expression = '"per_Envi_H" >75 AND "per_Envi_H" <100'
    else:
        Expression = "0"
    print Expression
    arcpy.SetParameter(2, Expression)

elif Environment_Index == "Percentage of Flooded Area":
    if Measure == "0-25%":
        Expression = '"Per_Flood" > 0 AND "Per_Flood" <25'
    elif Measure == "25-50%":
        Expression = '"Per_Flood" > 25 AND "Per_Flood" <50'
    elif Measure == "50-75%":
        Expression = '"Per_Flood" > 50 AND "Per_Flood" <75'
    elif Measure == "75-100%":
        Expression = '"Per_Flood" > 75 AND "Per_Flood" <100'
    else:
        Expression = "0"
    print Expression
    arcpy.SetParameter(2, Expression)
else:
    if Measure == "0-50kgs/yr":
        Expression = '"pest_kg" >0 AND "pest_kg" <50'
    elif Measure == "50-100kgs/yr":
        Expression = '"pest_kg" >50 AND "pest_kg" <100'
    elif Measure == "100-250kgs/yr":
        Expression = '"pest_kg" >100 AND "pest_kg" <250'
    elif Measure == "250-500kgs/yr":
        Expression = '"pest_kg" >250 AND "pest_kg" <500'
    else:
        Expression = "0"
    print Expression
    arcpy.SetParameter(2, Expression)

I am trying to get a string as an output parameter from the script tool and use the string as the expression for the select by attribute tool.

My user wanted to have a range of values to be selected from a field, I thought it was an easy process by writing a python code.

I tried using the set parameter for the output parameter. I don't understand where the problem is.
This is my script tool in the model builder


texts:
I am new to openlayers and I am using OpenLayers Version: 6.3.1.
I need to display a satellite map. I have used mapbox API to display the satellite map. It doesn't have streets and displays only the plain terrain. How do I create a new layer that displays the streets and country/region names with the satellite map?
I have attached my snippet for map below:

var raster = new TileLayer({
  source: new ol.source.XYZ({
    url: 'https://api.mapbox.com/v4/mapbox.satellite/{z}/{x}/{y}.png?access_token=My_ACCESS_TOKEN'
  })
});

var source = new VectorSource({wrapX: false});

var vector = new VectorLayer({
  source: source
});

var map = new Map({

  layers: [raster, vector],
  target: 'map',
  view: new View({
    center: [0,0],
    zoom: 0
  })
});


texts: Extracting, stacking and calculating Standard Deviation from NetCDF (.nc) to produce single standard deviation and Mean image in ArcGIS for Desktop?

texts:
I have a NetCDF file for several years of monthly Climate data ( Precipitation,PET,Temp,etc) data from CRU 3.22 (Climate Research Unit). I want to extract a subset of the dataset namely 2000-2010. And subsequently process the datasets into raster as annual mean and annual standard deviation and also produce a mean and a standard deviation for the period. To simplify the process i want to use either Python or ModelBuilder.

How to extract,stack and calculate Standard Deviation from a NetCDF (.nc) to produce a single standard deviation and Mean image?


texts: Apply feature class to feature class simulteaniously on multiple layers in Modelbuilder

texts:
I am trying to create a model in modelbuilder that takes in multiple layers (the amount determined by the person using the tool in the end) and as a first step these layers should by copied with the tool 'feature class to feature class' because the original layers are not allowed to be edited (they are from the company). However, feature class to feature class takes in only 1 input features.

Does anyone know how to perform this action on multiple layers at the same time?


texts: Is it normal that the values of a raster will change after reprojecting it?

texts:
I have a chlorophyll raster with WGS 1984 coordinate system and I reprojected it to UTM Zone 51 N using Project Raster tool and the nearest-neighbor interpolation option in ArcGIS 10.1.

The reprojected raster values changed (WGS1984: 0-58.7058 while the UTM Zone 51 N (reproject): 0-40.9507) after the procedure (see image below). What should I do? Is that really will happen after? The chl_reproject layer is the reprojected raster (UTM Zone 51N) from the original having a coordinate system of WGS 1984 (shown beneath chl_reproject).


texts: Why does Append on large annotation layer give WARNING 000594: Input feature falls outside of output geometry domains?

texts:
I am trying to Append about 3,000,000 records in an annotation layer to another, empty annotation layer. I'm using ArcCatalog and I get the following error messages:

WARNING 000594: Input feature 4657048: falls outside of output geometry domains.
WARNING 000595: c:\temp\annotation.fid contains the full list of features not able to be copied.

Approximately 100,000 features do not transfer. I have:

set the processing extent to Default and also Union of Inputs when running the Append tool
tried running this tool in a Python script, all the records will transfer but most of the field values do not.

At least some of the records that do not transfer are spatially in the midst of records that do transfer. Does anyone have some suggestions on how I can get all of the annotation records to transfer properly? Thanks & regards.

---- Update -----

When using the Append tool I set the schema type to test. Feel free to send along ideas, and I am going to try append the records in several pieces to see if that works. Thanks!

How did you create your empty annotation feature class?

I right clicked on my SDE database > new > feature class > annotation feature class

In particular what Coordinate System did you apply and how?

I applied the same coordinate system (State Plane Coordinate System, CA Zone V) as the annotation feature class I want to Append using the Import tool when I create the feature class.


texts: Merge DEM and Hillshade to create one file?

texts: I am wondering if I could merge DEM and Hillshade? What are the pros and cons for that?


texts: Is there a way to overlay a remote kmz file in google maps

texts:
I have a kmz file and want to overlay it on google maps. What I'm searching for is sth. like: https://www.google.de/maps?kmz=server/tour.kmz

I've already checked the query params for google maps, but couldn't find one for this purpose.

I don't want to upload and connect that kmz file with my google account.


texts: Displaying geographic information from Geoserver in Drupal through openlayers module

texts:
I'm trying to display geographic information coming from a local Geoserver instalation (local Tomcat server) in an also local Drupal installation (local Apache server) through openlayers module, but I think I am missing something.

Drupal is working on a local Apache server on port 8888; Geoserver is working on local Tomcat server on port 8080.

Consider the layer "Manhattan (NY) roads" which is installed as example in Geoserver. I can preview it in Geoserver (clicking on "openlayers" button) and the GeoJSON feed is working, if I select it from the dropdown menu, also in Geoserver's Layer Preview page.

GeoJSON url is:
http://localhost:8080/geoserver/tiger/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=tiger:tiger_roads&maxFeatures=50&outputFormat=application%2Fjson

and OpenLayers url is:
http://localhost:8080/geoserver/tiger/wms?service=WMS&version=1.1.0&request=GetMap&layers=tiger:tiger_roads&styles=&bbox=-74.02722,40.684221,-73.907005,40.878178&width=476&height=768&srs=EPSG:4326&format=application/openlayers

I tried using these URL's / portions of it in Drupal through the openlayers modules – creating a source, a layer and a map exposed through a block – but when I enable the block I see no map, only blank space.

I may be entering wrong urls, or selecting wrong source / layer types.

Can someone point me in the right direction?

EDIT 1

This is the Geoserver log, after trying to preview the map in Drupal:

2017-02-03 11:15:22,309 INFO [geoserver.wfs] -
Request: getServiceInfo
2017-02-03 11:15:22,314 INFO [geoserver.wfs] -
Request: getFeature
service = WFS
version = 1.0.0
baseUrl = http://localhost:8080/geoserver/
query[0]:
typeName[0] = {http://www.census.gov}tiger_roads
maxFeatures = 50
outputFormat = application/json
resultType = results
2017-02-03 11:15:22,315 INFO [wfs.json] - about to encode JSON


texts: How to convert to Australian Height Datum from Google Earth elevation?

texts:
I am trying to figure out how to convert from Google Earth elevation to Australian Height Datum (AHD).

Google Earth's elevation / altitude is based on
The vertical component (altitude) is measured from the vertical datum, which is the WGS84 EGM96 Geoid.
So I assume that WGS84 is Ellipsoidal value.
http://productforums.google.com/forum/#!topic/earth/fBar9afvRV0

The absolute difference between GDA94 and ITRF is now approximately one metre. Since the WGS84 reference frame is now aligned to the ITRF at the centimetre level, the absolute difference between WGS84 and GDA94 is also approximately one metre.
So I ignore the 1 metre difference.
http://www.dse.vic.gov.au/property-titles-and-maps/geodesy/geocentric-datum-of-australia-gda

And I use AUSGeoid09 http://www.ga.gov.au/ausgeoid/nvalcomp.jsp to compute.

I enter the latitude and longitude from Google Earth and into GDA94 Latitude and GDA94 Longitude text boxes. And I select Ellipsoidal for GDA94 Height (m) and enter the elevation value I get from Google Earth and click on compute. Then the next page will show me the AHD value for this location and assuming that the Google Earth elevation value I entered (WGS84) can have 1 metre difference (I don't know if this is + or - but I can try both +1 and -1).

I don't know if I am doing make any sense at all and would appreciate some helps.

And what I am trying to do is to compare the Q100 flood level with the AHD level on any given Australian point. Say the council give me a Q100 flood level in AHD, say 10m AHD, and I am trying to find out if the given point has a higher AHD value of 10 or not.


texts: Map projection classification - Tissot's indicatrices are perfect circles, but the projection is not conformal

texts:
Recently, led by the perfectly circular shape of the Tissot's indicatrices all over the projection plane, I have concluded that a particular projection is conformal. However, after calculating angular distortion, I have realized that it is actually not conformal.

How should I classify these projections, whose Tissot's indicatrices are perfect circles but the angular distortion is not zero (omega <> 0)?

I didn't name the projections above since none of them have yet been published, so their names are not meaningful yet. The first paper should probably be published next month, so I need to reassure that certain statements are correct before final publishing.

One of the projections is a spherical cube map projection used in Outerra engine, provided through the courtesy of Brano Kemen. The distribution of the texel aspect distortion can be seen in the paper Ellipsoidal Clipmaps. The texel aspect distortion is defined as the texel (a pixel in the texture space) width to height ratio after unprojecting to the surface of a planet. The distribution of texel aspect distortion and how the front face of the cube (actually a projection plane) looks like for the OSC is shown in the subfigure c) on the figures Err1 and SC10, respectively. I planned to give more links, but my reputation on this site (as a beginner) prevents me from doing so.

For everything I mention here, I'm using my own implementation. It is more flexible and for, at least, the order of magnitude faster than using any standard library (like Proj.4). Also, as I'm experimenting with the projections, using a standard library can be a limiting factor.

In the posted figures, I used the inverse projection of the elementary cells in the projected image and reprojected them to the surface of the Earth. After that, I measured the distance on the surface of the Earth along X and Y axis in the projection space. Here is the code that produces the aspect distortion metric:

double phi1, theta1, phi2, theta2;
inverse(x - delta, y, 0, phi1, theta1, false);
inverse(x + delta, y, 0, phi2, theta2, false);
double dLon = SphericalDistance2(phi1, theta1, phi2, theta2);

inverse(x, y - delta, 0, phi1, theta1, false);
inverse(x, y + delta, 0, phi2, theta2, false);
double dLat = SphericalDistance2(phi1, theta1, phi2, theta2);

double aspectDist = dLon / dLat;

The "inverse" function transforms projected (x,y) coordinates back to the surface of the planet (phi, theta). The last parameter defined whether an ellipsoidal model should be used. In the previous code, I'm using a spherical model, so the distance is measured on the spherical surface.

The problem is that the texel aspect distortion does not exist for the OSC but the angular distortion is not zero. I need to know how to reference such behavior of the projection. I thought it was conformity, but now I see that it is not. However, that property is very useful since causes no anisotropy in mapping projected image back to the surface of the planet. Or I made a terrible mistake.

For the angular distortion I'm using formulae from the Snyder’s "Map Projections – A Working Manuel" (pg.24, eq.4-10 through 4-14). Instead of partial derivatives, I'm using finite (symmetrical) differences in extremely small distances (~1e-10).