texts:
I'm just having a problem compiling postgis on a second machine, I wonder if anyone has any pointers:

./autogen.sh
./configure

gives

  PostGIS is now configured for i686-pc-linux-gnu

 -------------- Compiler Info -------------
  C compiler:           gcc -g -O2
  C++ compiler:         g++ -g -O2

 -------------- Dependencies --------------
  GEOS config:          /usr/local/bin/geos-config
  GEOS version:         3.3.0
  PostgreSQL config:    /usr/bin/pg_config
  PostgreSQL version:   PostgreSQL 8.4beta2
  PROJ4 version:        47
  Libxml2 config:       /usr/bin/xml2-config
  Libxml2 version:      2.6.26
  PostGIS debug level:  0

and make gives:

make -C liblwgeom
make[1]: Entering directory `/root/postgis-1.5.3/liblwgeom'
make[1]: Nothing to be done for `all'.
make[1]: Leaving directory `/root/postgis-1.5.3/liblwgeom'
make -C postgis
make[1]: Entering directory `/root/postgis-1.5.3/postgis'
Makefile.pgxs:17: warning: overriding commands for target `install'
/usr/lib/pgsql/pgxs/src/makefiles/pgxs.mk:92: warning: ignoring old commands for target `install'
Makefile.pgxs:63: warning: overriding commands for target `installdirs'
/usr/lib/pgsql/pgxs/src/makefiles/pgxs.mk:140: warning: ignoring old commands for target `installdirs'
Makefile.pgxs:82: warning: overriding commands for target `uninstall'
/usr/lib/pgsql/pgxs/src/makefiles/pgxs.mk:164: warning: ignoring old commands for target `uninstall'
gcc -O2 -g -pipe -Wall -Wp,-D_FORTIFY_SOURCE=2 -fexceptions -fstack-protector --param=ssp-buffer-size=4 -m32 -march=i386 -mtune=generic -fasynchronous-unwind-tables -I/usr/include/et -Wall -Wmissing-prototypes -Wpointer-arith -Wdeclaration-after-statement -Wendif-labels -fno-strict-aliasing -fwrapv -g -fpic -I/usr/local/include  -I/usr/include/libxml2 -I../liblwgeom -I. -I/usr/include/pgsql/server -I/usr/include/pgsql/internal -I/usr/include/et -D_GNU_SOURCE -I/usr/include/libxml2  -I/usr/include  -c -o lwgeom_accum.o lwgeom_accum.c
lwgeom_accum.c: In function ‘pgis_geometry_accum_transfn’:
lwgeom_accum.c:114: error: ‘WindowAggState’ has no member named ‘aggcontext’
lwgeom_accum.c:308:66: error: macro "InitFunctionCallInfoData" passed 6 arguments, but takes just 5
lwgeom_accum.c: In function ‘PGISDirectFunctionCall1’:
lwgeom_accum.c:308: error: ‘InitFunctionCallInfoData’ undeclared (first use in this function)
lwgeom_accum.c:308: error: (Each undeclared identifier is reported only once
lwgeom_accum.c:308: error: for each function it appears in.)
make[1]: *** [lwgeom_accum.o] Error 1
make[1]: Leaving directory `/root/postgis-1.5.3/postgis'
make: *** [postgis] Error 2


texts:
I have a set of postcode boundary data that currently goes down to the most detailed level possible, like this:

AB1 2CD

I would like to 'simplify' the data so that it just shows the first part of the postcode, and thus a simpler, higher-level boundary. Essentially I just want:

AB1

In MapPublisher, I can import the data, apply a Selection Filter that looks for 'AB1' (for example), and then use Illustrator's pathfinder tool to merge this into one polygon consisting of just the outer boundaries of the smaller constituent polygons.

However, the problem I have is that MapPublisher can't export this data out as a Shapefile.

In QGIS, I can select individual sectors by simply searching for 'AB1' in the attribute table. However, I don't know how to go about A)'merging' the individual areas into the one big shape or B)exporting them all to a new Shapefile.


texts: Is there cartographic convention for displaying text label of highest and lowest values in legend?

texts:
The first I have found is to have the lowest value and highest value be displayed like the legend shown below (image source)

The second is to have relative numbers at the top and bottom (image source)

Should I display 1-100 items or less than 100 items? Is one method better than or preferred over the other?


texts:
I am trying to convert an .asc file to .dwg using QGIS but I receive the following error:

  ERROR 1: CAD driver has no creation capabilities

Could anyone please advise what the problem could be?

I suspect it is a bug.


texts: raster doesn't redraw to fit the canvas after a zoom in and a pan

texts:
i am working with qgis 2.8.2, c++ api and qt 4.8.5
i have a layer in the mapcanvas, this layer contains a raster ( a georef tif picture )
when i zoom in, the raster fill the entire canvas , ok but if i grad it in any direction, in the left bottom corner by exemple, the rest of the canvas stay empty
i want to know how i can redraw the canvas in order that the zoomed raster fit all the canvas and a little region on the screen

i try to use updateMap(); in the mouseReleaseEvent but nothing change


texts: Error installing Rasterio in PyCharm and in CMD Line

texts:
Error installing Rasterio in PyCharm (Settings > Project Interpreter)

Collecting rasterio
  Using cached https://files.pythonhosted.org/packages/1a/d5/6cad936ff9beb89e79b1dae946f3097c70f7ab71cda5a9e356aec2e9ce9b/rasterio-1.0.8.tar.gz
    Complete output from command python setup.py egg_info:
    INFO:root:Building on Windows requires extra options to setup.py to locate needed GDAL files. More information is available in the README.
    ERROR: A GDAL API version must be specified. Provide a path to gdal-config using a GDAL_CONFIG environment variable or use a GDAL_VERSION environment variable.

    ----------------------------------------

Command "python setup.py egg_info" failed with error code 1 in C:\Users\X\AppData\Local\Temp\pycharm-packaging\rasterio\
You are using pip version 9.0.1, however version 18.1 is available.
You should consider upgrading via the 'python -m pip install --upgrade pip' command.

Same error when trying to install in command line; following is cmd line entry:

cd C:\Python27\ArcGIS10.6\Scripts
pip freeze
pip install rasterio

Here's the command line return:

Collecting rasterio
  Using cached https://files.pythonhosted.org/packages/1a/d5/6cad936ff9beb89e79b1dae946f3097c70f7ab71cda5a9e356aec2e9ce9b/rasterio-1.0.8.tar.gz
    Complete output from command python setup.py egg_info:
    INFO:root:Building on Windows requires extra options to setup.py to locate needed GDAL files. More information is available in the README.
    ERROR: A GDAL API version must be specified. Provide a path to gdal-config using a GDAL_CONFIG environment variable or use a GDAL_VERSION environment variable.

    ----------------------------------------
Command "python setup.py egg_info" failed with error code 1 in c:\users\starb\appdata\local\temp\pip-install-dw_0ma\rasterio\

I can't find the gdal-config path. I guessed gdal-data and added that to Windows Environments, but no dice. I installed GDAL following these instructions: https://sandbox.idre.ucla.edu/sandbox/tutorials/installing-gdal-for-windows


texts: qgis-mapserver: GetCapabilities returns no layer

texts:
My problem is similar to
this thread

I've installed QGIS 2.4 and QGIS server 2.4 on Fedora. I created a simple project containing a vector layer (shapefiles)
GetCapabilities request works but there is no layer in it. Except the Layer tag containing the name of the project.

Example:

<Layer queryable="1">
<Name>test3</Name>
<Title>test3</Title>
<EX_GeographicBoundingBox>
<westBoundLongitude>-68.5</westBoundLongitude>
<eastBoundLongitude>-68.5</eastBoundLongitude>
<southBoundLatitude>44</southBoundLatitude>
<northBoundLatitude>44</northBoundLatitude>
</EX_GeographicBoundingBox>
<BoundingBox CRS="EPSG:32198" maxx="0" minx="0" maxy="0" miny="0"/>
</Layer>

here is an example of my query

http://localhost/qgis/qgis_mapserv.fcgi?map=test3.qgs&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities

I tried in the browser and curl command line and got the same result.

All the characters in the path are ascii.
Perhaps it is a privilege access problem but I could not resolve it. My path is readable by all, since GetCapablities works it can read my .qgs.

Any idea on what the problem is?


texts: Customising line thickness of each line mapped in R using tmap package?

texts:
I have generated a map showing the intensity of burglaries that have occurred on a street segment. I have used the tmap package from R to generate the desired street-level map. The codes are as follows:

   map.burglary =  tm_shape(street.crime.rates.lines) +
                   tm_lines(col = "b_rates",
                            style = "fixed",
                            breaks = c(0, 1, 10, 20, 50, 500),
                            labels = c("<1.0", "1.0-9.9", "10.0-19.9","20.0", 49.9", "50.0+"),
                            textNA = "Unavailable",
                            colorNA = "gray87",
                            title.col = c("Burglary (per 10 households)")) +
                   tm_compass(type = "8star", position = c("right","top")) +
                   tm_scale_bar(breaks = c(0, 0.25, 0.5, 0.75, 1)) +
                   tm_layout(frame = TRUE, legend.title.size = 0.8, legend.outside = TRUE, legend.outside.position = "right")
    map.burglary

The resulting output:

I want to increase the line thickness of each line as the intensity of burglaries represented on the map becomes larger. I want this detail to be reflected in the map and legends. Does anyone know how to code this in R using the tmap package?


texts: Non-alignment of vector and raster layers in QGIS

texts:
I have a vector map layer, raster climate data and added a delimited text points layer. My problem is that the raster extent does not quite match the edge of the vector map, and as a result, some of my coastal data points are being excluded when I sample raster values.
I've added a screen shot (hopefully it will be clear from this). I have made sure that the CRS projections are the same. Is there anything else I can do?


texts: How can my client catch mapservice timeouts?

texts:
I've got a gp service that exports pdfs using connections to mapservices.  I'd like to write to a log file when a request to a service times out.  For WMS mapservices I see I can use IWmsClientEvents outbound interface on the WmsConnection, but for the AGSServerConnection I don't see any outbound events.

Is there a way to find out when a draw request to an arcgis server mapservice times out?


texts: How send parameters to a Dockable Window in ArcObjects?

texts:
I need send a IFeature element to a dockable window.

I tried creating a public static variable, in eManager Class, that is inicialized with a feature that catch a function:

public static IFeature pFeature;
    internal void ShowEntityDetails(IFeature feature)
    {
        pFeature = feature;
        IDockableWindow dockWindow = GetDockWindow(ThisAddIn.IDs.GMREAddin_Modules_Model_View_vEntityDetails);
        dockWindow.Show(true);
    }

And I get it in dockable window construct of this way:

public vEntityDetails()
    {
        InitializeComponent();
        IFeature pFeature = eManager.pFeature;
        if (pFeature == null) return;

        string streetName = pFeature.Value[pFeature.Fields.FindField("ALIAS")].ToString();
        MessageBox.Show(streetName);
    }

It only works the first time.

There any more professional way to do that?


texts: I want to select a polygon using QGIS. The polygon is only represented via it's boundaries (I chose not to fill the area in order to highlight the problem). Why are not all lines that surround the polygon selected? See picture below. A similar problem occurs when two polygons overlap: only one polygon is highlighted.


texts: Using centroid point to generate polygon using ArcPy?

texts:
Is this possible to generate polygons from points (which are centroids of future polygon)? I got shapefile with ~250 points + they have got two attributes: length and width.

I don't need complete workaround I need tip on how this can be achieved.

Below is sketch picture of what I want to achieve.


texts: Replacing action triggered script by one supplied through QGIS plugin?

texts:
In ArcView 3 all controls (menu choices, buttons, tools) had scripts associated with them that could be changed simply by supplying a script with the same name through an extension (plugin). That way you spared designing a new control and could keep your panels tidy while improving and extending functionality of the program.

Now I'd like to do the same thing with QGIS: E.g. I'd prefer to preselect the default directory of the "add ogr layer" action by the path of the first active layer in my map.
So if it's possible I'd like to replace or modify the mActionAddOgrLayer.triggered Script, as long as my plugin is loaded.

Another way to achieve this would be to modify the built-in system response to a change of the active layer(s) by changing the variable that by default saves the directory of the last file added. I'm sure there is a signaling but I don't know if and how responses can be changed or added without adding additional gui-controls (actions).


texts:
I have ArcGIS 10.2 for Desktop and I have created code which does the following. Abstracts from 1st row 1st column 10, (1500-10) adds it at the second and then repeats A-=A where A= 1500. I want to write a tool with that script that will ask from user to put the fclass and then do the procession but I have no idea how to do that altough I have read several things about creating tools. Propably I must use get parameter as Text but I really don't know how.

Q      A
1500   1500
1340   1490
1234   1480
1543   1470

import arcpy
f = "c:/Wi/Arm/Fo.shp"
field = ["Bo", "Cond"]
i = 10
with arcpy.da.UpdateCursor(f, field, "FID = 0") as cursor:
    for row in cursor:
        row[0] = row[0] - 10
        cursor.updateRow(row)
        Bo = row[0]
        Cond = Bo
        cursor.updateRow([Bo, Cond])
with arcpy.da.UpdateCursor(f, ['Cond'], "FID >= 1") as cursor:
    for row in cursor:
            Cond = (Cond - i)
            cursor.updateRow([Cond])


texts:
We are working on a project where we want to make a route planner that takes into account the amount of sun on the route. The goal would be to calculate 'solar car' friendly routes. We are using PostGIS, because of the pgRouting extension for the route finding.

We have building polygons with a height attribute available. So far we have not been able to find PostGIS functions that could deal with calculating shadows.

How do you calculate polygon shadows in PostGIS?

Using PostGIS to generate building shades seem to be simply wrong (correct us if we're missing anything). The answer suggests extruding a 3D polyhedral from a 2D surface. That does not seem to generate a shadow. Only if you visualize it in 2D.

This is what we want as result (picture from H2GIS):


texts: Plotting gam spatial estimates over map?

texts:
I have imported Italy map with the following code:

library(ggplot2)
library(mgcv)
library(ggmap)
Country <- "Italy"

Get_Map_Country <-
  get_map(
      location = Country
    , zoom = 6
    , scale = "auto"
    , maptype = "satellite"
    , messaging = FALSE
    , urlonly = FALSE
    , filename = "ggmapTemp"
    , crop = TRUE
    , color = "color"
    , source = "google"
    , api_key
    )

Country1 <-
  ggmap(
      ggmap = Get_Map_Country
    , extent = "panel"
  #  , base_layer
    , maprange = FALSE
    , legend = "right"
    , padding = 0.02
    , darken = c(0, "black")
    )

Furthermore, I fitted a gam model in which the response is assumed to be Poisson
and the predictors are longitude and latitude.

mod1 <- gam(Value ~ s(Longitudine, Latitudine), data = data2017, family = "poisson")

Now, how can I overlay the results to the ggmap above?
For instance, I would like to obtain something like:

obtained by using the "mfasiolo/mgcViz" package installed with devtools from github

library(devtools)
install_github("mfasiolo/mgcViz")
library(mgcViz)

b <- getViz(mod1)

pl <- plot(sm(b, 1)) + l_fitRaster() + l_fitContour()
pl

R base graphics analogous version can be obtained by:

plot(mod1, scheme = 2)

But good-looking (for publication), slightly transparent over the italian map.

How can I obtain this?

Simple reproducible example

> y <- c(10,   1,   2,   2,   2,   0,  19,   3,  2,   1,   7,  16,   5,
> 13,  81,   1,   9,  96,   2,   9,  3,   4,   0, 1,  12)
>
> long <- seq(9, 13, length.out = 25)
>
> lat <- seq(38, 43, length.out = 25)
>
> mod1 <- gam(y ~ s(long, lat, k = 3), family = "poisson")
#Now, I need not overlay Country1 and something like plot(mod1)


texts: Comparing field definitions in Access DB to GDB using ArcPy?

texts:
New to the GIS world, but have a little background with Python. I'm attempting to compare an Access DB to a GDB and update the GDB for specific fields if there is a difference. I've read multiple places that updating a field like length or width directly with SQL can cause issues, so I wanted to make sure I'm using the appropriate modules to update the fields. I found a similar link below that is attempting to something similar, but I'm only wanting to update the differences in specific fields. The other issue is that the field columns are named differently, utilizing SQL I could just rename the column and match the values appropriately, but I haven't found any documentation stating that I can do the same with modules such as TableCompare_management, and join.

Just to confirm, can I not just compare each row from the access db by a unique ID to the GDB, and update the fields that are different? Or will that really break things?


texts: Making a screenshot from QGIS2web Openlayers map

texts:
I would like to make a screenshot from my Openlayers map
My code so far looks like this:

The map comes from qgis2web.js  file.
The console says nothing, but the map is completely blank when clicked at: generate screenshot button.
Why can't I grab the Qgis2web map content in  section properly?


texts: Creating quintiles in QGIS when the number of observations per row is variable

texts:
I have a list of 100s of schools (1 per row), and a variable column with the total students in each school. The number of total students in each school varies widely.
I want to split the list of schools into five categories such that each category contains roughly the same number of STUDENTS; the number of schools in each group doesn't matter.

Can I do this in QGIS?


texts: MapServer Cluster, problem with GetFeatureInfo

texts:
I'm new with MapServer and I tried to use as cluster.
When I click on the map to see info about my points, nothing pop up. I can see info about different layers but not for the one with the cluster.
Layer with the cluster:
LAYER
NAME "eqp_point_collecte"
GROUP "eqp_dechet"
TYPE POINT
TEMPLATE "ttt"
CONNECTIONTYPE postgis
PROCESSING "CLOSE_CONNECTION=DEFER"
CONNECTION "user=${PGUSER} password=${PGPASSWORD} dbname=${PGDATABASE} host=${PGHOST} port=${PGPORT}"
DATA "geom FROM (select * from ${PGSCHEMA_GEODATA}.eqp_point_collecte) as foo USING UNIQUE gid USING srid=2056"
METADATA
    "wfs_title" "eqp_point_collecte" ##REQUIRED
    "wfs_srs" "EPSG:2056" ## REQUIRED
    "wfs_enable_request" "*"
    "wfs_include_items" "all"
    "wms_title" "eqp_point_collecte" ##REQUIRED
    "wms_srs" "EPSG:2056" ## REQUIRED
    "wms_enable_request" "*"
    "wms_include_items" "all"
    "gml_include_items" "all"
    "ows_enable_request" "*"
    "ows_geom_type" "point" # For returning geometries in GetFeatureInfo
    "ows_geometries" "geom" # For returning geometries in GetFeatureInfo
    "ows_title" "eqp_point_collecte"
    "gml_types" "auto"
    "gml_featureid" "gid" ## REQUIRED
    "gml_geometries" "geom"
    "gml_geom_type" "point"
END
DUMP TRUE # For GetFeatureInfo
STATUS ON
PROJECTION
  "init=epsg:2056"
END
CLASSITEM "Cluster_FeatureCount"
EXTENT 2470000 1110000 2840000 1200000
#SYMBOLSCALE 500
MAXSCALEDENOM 25000
###########################
# Define the cluster object
###########################
PROCESSING "CLUSTER_GET_ALL_SHAPES=ON"
PROCESSING "ITEMS=gid"
CLUSTER
  MAXDISTANCE 30
  ##########################################################
  ### TIPS : MAXDISTANCE > SIZE du symbole évite les bugs###
  ##########################################################
  REGION "ellipse"
END
################################
# Class1: For the cluster symbol
################################
CLASS
    NAME "Déchetterie"
    EXPRESSION ("[Cluster_FeatureCount]" != "1" AND '[type_collecte]' = 'déchetterie')
    STYLE
       SIZE 20
       SYMBOL "déchetterie"
       COLOR 255 0 0
    END
END
CLASS
    NAME "Eco-point"
    EXPRESSION ("[Cluster_FeatureCount]" != "1" AND '[type_collecte]' = 'écopoint')
    STYLE
       SIZE 20
       SYMBOL "éco point"
       COLOR 255 0 0
    END
END
################################
# Classes: For single stations
################################
CLASS
    NAME "PET"
    EXPRESSION ('[type_dechet]' = 'PET')
    STYLE
       SYMBOL "PET"
       SIZE 20
       OUTLINEWIDTH 0.05
       OUTLINECOLOR 240 175 55
       COLOR 255 255 255
    END
END


texts:
I'm writing this code to write the NDVI value and the date where picture was taken into my attribute table, later i'll export to excel. The code is downloading the images with no problem, but when it pass through the raster layer, it retrieves me an error! I don't understand what is going wrong...

import arcpy
import os
import datetime
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
NDVI_tempo = 'Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL_PL'
arcpy.env.workspace = 'Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\NAL.gdb'
arcpy.SplitByAttributes_analysis("Area_produtiva3", NDVI_tempo, ['ID_TALHAO'])
#Inicio do download das imagens

down = 'Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL_PL'
image_server = 'Y:\\E - IMAGEM SERVER\\IMAGEM_SERVER.gdb\\F01_PLANETLAB'

arcpy.env.workspace = r'Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL_PL'
fcs = arcpy.ListFeatureClasses()

for fc in fcs:
    arcpy.DownloadRasters_management (image_server, down, "", fc, "CLIPPING", "CONVERT_AS_REQUIRED", "TIFF", "NONE", "", "NO_MAINTAIN_FOLDER")

    rasters = arcpy.ListRasters()
    for tiff in rasters:
            arcpy.MakeRasterLayer_management(tiff, "Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\temp_rast1.tif", "", "", "3")
            arcpy.MakeRasterLayer_management(tiff, "Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\temp_rast2.tif", "", "", "4")

            arcpy.CopyRaster_management("Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\temp_rast1.tif", "copy_rast1", "","", "", "", "", "", "", "", "", "")
            arcpy.CopyRaster_management("Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\temp_rast2.tif", "copy_rast2", "","", "", "", "", "", "", "", "", "")

            Property1Result = arcpy.GetRasterProperties_management ("copy_rast1", "MINIMUM")
            Property2Result = arcpy.GetRasterProperties_management ("copy_rast2", "MINIMUM")

            Property1 = Property1Result.getOutput
            Property2 = Property2Result.getOutput

            Correcao1 = Minus ("copy_rast1", "Property1")
            Correcao2 = Minus ("copy_rast2", "Property2")

            Rastdiv1 = Minus ("Correcao2", "Correcao1")
            Rastdiv2 = Plus  ("Correcao2", "Correcao1")

            NDVI_tiff = Divide (Rastdiv1, Rastdiv2)
            Zonaltif = arcpy.ZonalStatistics (fc, ['ID_TALHAO'], NDVI_tiff, "MEAN", "NODATA")
            Zonalper = Times (Zonaltif, 100)
            Zonalint = Int (Zonalper)
            arcpy.RasterToPolygon_conversion (Zonalint, Zonalpoly, "SIMPLIFY", "VALUE")
            arcpy.Intersect_analysis ([Zonalpoly, fc], "Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL_PL\\teste\\Intertif", "ALL", "", "")
            arcpy.AddField_management (Intertif, Data, "TEXT", "", "", "", "", "", "", "")
            tiff_name = arcpy.Describe(tiff).name
            field_name = "{}".format(tiff_name[5:13])
            field_value = datetime.datetime.strptime('field_value', '%Y%m%d').date()
            arcpy.CalculateField_management (Intertif, Data, field_value, "", "")
            arcpy.Delete_management (tiff, "")

Traceback (most recent call last):
  File "Z:\04 - Temp\Charles\NDVI_tempo\NAL\NAL.gdb\NDVI_RF#Script3_NDVIRF.py", line 40, in
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\management.py", line 6219, in MakeRasterLayer
    raise e
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\management.py", line 6216, in MakeRasterLayer
    retval = convertArcObjectToPythonObject(gp.MakeRasterLayer_management(*gp_fixargs((in_raster, out_rasterlayer, where_clause, envelope, band_index), True)))
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing_base.py", line 506, in
    return lambda *args: val(*gp_fixargs(args, True))
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
 ERROR 000865: Input raster: 20170312_124653_0e1f_381_s65.tif does not exist.
Failed to execute (MakeRasterLayer).
Failed to execute (Script3).


texts: Choosing Projection for Gulf coast region of USA?

texts: What are the most appropriate datum and projection I could use for the Gulf coast region - spanning  Louisiana, Mississippi, Alabama, and Florida? I have to analyze block group data along the coast.


texts: Is there a standard as to how many meters or feet if a GPS number (latitude and longitude) is off by a tiny amount, such as 0.0001?

texts:
For example, if the latitude and longitude is

25.1054, 121.5973

then what if it is off by a small amount, such as 0.0001:

25.1055, 121.5974

Then can it readily tell how many feet or meters it is off by?


texts:
I am trying to follow the http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/ but it seems to be out of date, for example many of the raster methods no longer work as indicated or even exist (eg drawingStyle).
Can anyone suggest a better online resource for learning this stuff?


texts: How do I import OS Sites Layer into QGIS/PostGIS?

texts:
I've found the OS Translator II plugin which works perfectly for the Topology Layers, but it doesn't have an option to import OS Sites Layers.

How can I import this product into QGIS/PostGIS? I've tried naively importing the gml.gz file but it only gives the following fields.

id
geom
gml_id
toid
version
versionDate
reasonForChange
refToFunctionalSite
accessType
refToRoutingPoint
natureOfAccess
refToITNRoadNode
itnRoadNodeVersionDate

Answers using the QGIS GUI or Python psycopg2 are equally useful.


texts: Custome QGIS application with WPS and GDAL/OGR

texts:
I have basic question on how to start a project with QGIS. There are two kinds of users, and they are supposed to sign in (Login in) before they start using so that depending on the user type the QGIS interface and plug-ins will be presented to them. So basically there are two different QGIS interfaces.
So my questions are:

Is there any way I can add a login interface before QGIS starts?
Can we customize the QGIS in different ways so that different users can have access to different QGIS interface and plug-in?

After login the first user is supposed to do editing, and can save it on a local disk to continue later. But after some certain button is clicked, the process continue to the second user whose main job is to approve the edited features done by the second user. If approved, the newly edited data is committed to a database. If not the it will be reverted back to the first user with some comments.

Is there any way that such a process can be implemented in QGIS? Would WPS be any help?
If there is a need in the process to change one format to another, could we use the GDAL/OGR library?

I heard and most examples that I come across show how to convert formats with GDAL/OGR using command line.

Is there any possible way you could use the API in your code? Like with in your WPS code or inside QGIS?

I am guessing some of my questions may not give sense since I am all new to open source based GIS programming, but I would really be grateful for any kind of help.

You be cool!


texts: Why modified feature resets it's position on zoom in/out in Open Layers 6

texts:
I have used select and modify interactions with OpenLayers 6, used feature option for modification. These works fine but face an issue with modify interaction. When I change the position of a feature using modify interaction it reflect updated position on map, but as I zoom in or out the modified feature resets it's position to original(which was before modification).
This fails when cluster layer is used, when no cluster layer it works fine. What might be the issue here.


texts: Sentinel-2 and Landsat-8 future acquisitions

texts:
In order to schedule in-situ measurements at the time of Sentinel-2 and Landsat-8 overpass for my AOI, I would like to access their future acquisitions plans for the summer 2020 (July/September).

I have read in this post that the list of all future acquisitions is difficult to predict. I would like to better understand why: if possible I would like to build a preliminary schedule just to have an idea, even if it's not 100% reliable.

So far I have the following sources of acquisition plans:

ESA publishes the Sentinel-2 acquisition plan, but its time coverage is only a couple of weeks
USG publishes a Landsat acquisition tool which seems to fit my needs (it displays overpasses up to 2022)

The Landsat Acquisition page explains that the acquisition plan is based upon seasonality, land definition, historical cloud cover, gain settings, and sun angle.

How come the Landsat-8 tool has such a long time coverage and the Sentinel-2 tool doesn't?
Should I not trust the Landsat-8 tool?
Are the acquisition plans of both satellites not subject to the same uncertainties?


texts:
I have connected my ArcMap to my PostgreSQL db using the Create Enterprise Geodatabase tool in ArcMap. I have multiple feature classes and the relevant relationships between them loaded into the db correctly.

However, I am attempting to add a unique sequence to the feature classes using the postgres sequence function as below:

CREATE SEQUENCE my_sequence;

I then set the default for an integer colum in my tables and assign the sequence as below:

SET DEFAULT nextval('my_sequence');

However, the sequencer works perfectly if I insert data inside PostgreSQL using a SQL statement. If I edit the same table in ArcMap I get the following error when trying to save my edits:

ERROR: duplicate key value violates unique constraint
DETAIL: Key(uniq_id)=(0)already exists

Where am I overlooking something?


texts: Earth Engine: downloaded image contains only 1 pixel after mosaic with python version

texts:
I followed the tutorial of getting a free-cloud landsat8 image from the official site with python version. But I got only 1 pixel with the function getDownloadUrl.Here's the code

def maskClouds(image):
    scored = ee.Algorithms.Landsat.simpleCloudScore(image)
    return image.updateMask(scored.select(['cloud']).lt(20))

def addQualityBands(image):
    return maskClouds(image).addBands(image.metadata('system:time_start'))

landsat = ee.ImageCollection('LANDSAT/LC08/C01/T1').filterDate(2017-08-01, 2018-08-01)
landsat = landsat.map(addQualityBands).qualityMosaic('system:time_start')
landsat1 = landsat.select('B1','B2','B3','B4','B5','B6','B7')
url = landsat1.getDownloadUrl({
    'region': [[-120, 35], [-119, 35], [-119, 34], [-120, 34]]
})


texts: Fixed Distance Buffer Gives Inaccurate Distances

texts: I have a world map shapefile showing all continents. Each of the shapefile is projected and I need to buffer each of the continents with a 1000 km buffer, but all my buffers show incorrect distances near the polar regions. They appear to be fixed-width as opposed to fixed-distance. They should look distorted near the poles. Am I missing something?


texts: QGIS 3.10: why is PCS project, loaded with a single PCS layer, showing only geographic (angular, lat/long) coordinates instead of linear ones

texts:
[OP: my error, circle CRS in LR corner is Project's CRS, not layer's.]
QGIS 3.10; Ubuntu 18.04
Glossary for new users: CRS-Coordinate reference system (one of two types, geographic or projected); GCS-geographic coord. system, on-the-globe, lat/long, angles; PCS-projected coordinate system, onto a 2D map, lengths
Here are my steps:

start QGIS and set Project->Properties->CRS to a PCS (EPSG:5070)
load a PCS raster file as a layer (saved from a previous GCS raster layer, as in this thread)
verify that the layer's CRS is identical to the project's
note that the coordinates displayed below the layer view are still GCS, lat/long, only.

Another thread advises to edit the attribute table, but that applies to vector layers only.
I know I've expressed a considerable amount of thickness regarding how to handle mixed layers of GCS and PCS coordinates over some years, but it's beginning to feel like I'm closing in on it.
Q. Is there a 'toggle' between displaying GCS and PCS in the main QGIS window?  No, but discussion follows.
Response to @Ian Turton's request for 'gdalinfo layer-file'
    grdn36w112_1$ gdalinfo w001001.adf
Driver: AIG/Arc/Info Binary Grid
Files: .
       ./sta.adf
       ./log
       ./w001001.adf
       ./hdr.adf
       ./prj.adf
       ./metadata.xml
       ./dblbnd.adf
       ./w001001x.adf
Size is 3612, 3612
Coordinate System is:
GEOGCS["NAD83",
    DATUM["North_American_Datum_1983",
        SPHEROID["GRS 1980",6378137,298.257222101,
            AUTHORITY["EPSG","7019"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6269"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4269"]]
Origin = (-112.001666666700004,36.001666666664136)
Pixel Size = (0.000277777777778,-0.000277777777778)
Corner Coordinates:
Upper Left  (-112.0016667,  36.0016667) (112d 0' 6.00"W, 36d 0' 6.00"N)
Lower Left  (-112.0016667,  34.9983333) (112d 0' 6.00"W, 34d59'54.00"N)
Upper Right (-110.9983333,  36.0016667) (110d59'54.00"W, 36d 0' 6.00"N)
Lower Right (-110.9983333,  34.9983333) (110d59'54.00"W, 34d59'54.00"N)
Center      (-111.5000000,  35.5000000) (111d30' 0.00"W, 35d30' 0.00"N)
Band 1 Block=256x4 Type=Float32, ColorInterp=Undefined
  Min=1079.405 Max=3847.728
  NoData Value=-3.4028234663852886e+38


texts: Splitting GeoTIFF into tiles in QGIS NextGIS

texts:
There is a GeoTIFF of around 3.5 Gb. I would like to get 6 sliced jpg tiles from it (width is one second of the original file, and hight is one third) with JGW bindings.
How is it possible to implement in QIGS?


texts: Creating mosaic dataset in Python without using ArcGIS for Desktop Standard/Advanced?

texts: Is there a way to create mosaic dataset via Python, numpy or gdal libraries without having a Standard/Advanced ArcGIS?


texts: How to get GPS location in OSMDroid without Internet connection?

texts:
I made a little app in Android using OSMDroid library, my map is loading from assets. Now I want to put my GPS location on that map, but without using Internet connection. I found these

GpsMyLocationProvider provider = new GpsMyLocationProvider(getApplicationContext());
        provider.addLocationSource(LocationManager.GPS_PROVIDER);
        MyLocationNewOverlay myLocationNewOverlay = new MyLocationNewOverlay(provider, mapView);
        myLocationNewOverlay.enableMyLocation();
        Toast.makeText(getApplicationContext(),"my location "+        provider.getLastKnownLocation()
                ,Toast.LENGTH_LONG).show();
        mapView.getOverlays().add(myLocationNewOverlay);

but nothing changes.


texts: Determine polygon mid point, align samlpling plots by cardinal points and distance to edge (Esri ArcGIS 10.2)

texts:
My sampling design requires me to do a couple of things:

First I need to determine the mid-point of an irregular polygon feature. Then I'd need to create a layer of 4 line features, that intersect in the mid-point of said polygon and are all oriented N-S, NE-SW, E-W and SE-NW. As a third and final step I need to highling a point on each line that lies exactly halfway between the midpoint of the polygon and the edge of it, plus one additional point at the very center. In the end, it should look something like this:

Is this feasible at all? Thing is, I have to do this for 32 polygons and it's just a small step for my project, so I want to spend as little time for it as necessary and thus need a method that works fast. The lines are actually just a tool to get to the 9 plot points, as plots determine the locations where I will do my sampling. So if anyone has any idea to get my 9 sampling-spots to be aligned the way I described, I'm open for all suggestions.


texts: How to get the layers of my project to show while tracking with GPS?

texts:
I've been using/testing several versions of various GIS software and run into occasional issues. Right now, my biggest with QGIS is with GPS connection/downloading from/tracking. These are all critical issues for me.

1.) I think I did finally figure out the need to set data transfer to NMEA in/out on my Garmin GPSMap60CSX. QGIS repeatedly crashed on clicking connect, I hope this was my only error and is now corrected.

2.) Downloading from/uploading to GPS unit... clicked to download tracks and waited 10 minutes before giving up and canceling with no success. Any tips/pointers?

3.) I've posted in the old forum about this one... I am unable to see any layers of the open project when live tracking. Last time (I think using Copiapo) I was eventually able to zoom until the .shp file was visible (basically stumbled onto that one), but not the .sid file that formed the base map. Now, using Wroclaw there is nothing at all visible except the gps track. Just creating a track doesn't help, I need to be able to see my location in relation to various physical features and property boundaries in real-time.


texts: Merge and count adjacent polygons in QGIS

texts:
I have a shapefile that shows single steps as individuals shapes and im trying to figure out the total count of adjacent steps. The little model I came up for now includes using the v.dissolve function from the GRASS-library to dissolve all shapefiles that are adjacent and then with join attributes by location to get the count of all shapefiles that touch, intersect or are inside the dissolved polygon. For most of the shapefiles this works fine, but as so often, it doesnt work fine on some others for different reasons.

Working fine:

Not working fine:

1: There may be some digitalization error in the first problematic example. The amount of steps I really want is 51 but Im getting 78, because the cross bar and the little rectangular shape inside the steps divides single steps into several steps which are then counted separately.

2: In the second example two stairs seem to be next to each other resulting in 49 individual stairs. But actually I would like to have 2 results, a step count for each stair. The left stair should have 27 steps and the right stair has only 22 steps.

3: The third example results in 66 steps. Maybe some steps are missing at the inside and actually there would only be 3 stairs next to each other but for now I would like to get 4 individual results.

4: In example 4 the wrong step count is 43. Two shapefiles are inside the stair on the left, dividing single steps further and are counted twice then. The correct step count for left stair would be 25.

So i have to adapt the model to take those exeptions into account but im not quite sure how to achieve this in the best way. I thought about building centroids of all single steps, joining only close points together as a path or by buffered area and then counting the steps that intersect with this path or that buffered area. This would maybe solve problems 2 & 3 but could be problematic with 1 & 4.

Another idea was to restrict single steps to have a minimum size but this would not affect problem 2 & 3.

And of course there would be the possibility to combine the above mentioned methods to solve all problems. But before starting to implement that method and see if it really works and what new problems might arise, I wanted to ask you guys for some experience or thoughts on this :)
Maybe there is some plugin that im not aware of or another process which could be usefull for this case.

Short update: The idea with the centroids, that are buffered, dissolved and then intersected with the steps to get their count seems to solve some of the problems (2 & 3 for example) but some problems remain. Of course..

And here are some examples, where it still doesnt work correctly.

The main problem is gonna be the exact parametrization of the buffer size. If the buffer size is bigger it would probably solve the brown stairs, be complicated with the blue one and could definitly mess up the stairs on the right, so that both slim stairs will be regarded as one, which might btw also be the case, as they are only about 0,3 meters broad. But it would then double the amount of steps in that buffered area.

The problem with the green buffered area, that two of such buffered areas are inside one stair should be solved quite easily.

Question: How can I calculate the max with and max length of every single shapefile, to assign this to a variable buffer size?


texts: Importing AUTOCAD .dwg to Arcgis Symbols

texts:
My client has some symbol drawn in autocad.. they are saved in .dwg files.

How can I import them into ArcGIS symbols?

I have FME 2013, I am not sure if that helps.


texts: Find gravity normal by ECEF or GPS coordinate of point

texts:
How to find gravity normal by ECEF or GPS coordinate of point?

Is it the same as normal to elipsoid?

Is this ordinary way to calculate normal to ellipse?

I'm not quite sure which altitude is GPS use, but it seems if it use geodetic latitude it already perpendicular to elipsoid.

And same about ECEF

To compute normal(I want it in ECEF) I need 2 points, 1st point I get directly from GPS coordinate converting it to ECEF and 2nd point I need to compute.

For example point that have Z=0.

code :

Vec4d ecef= GPStoECEF(gps);

double Radlat = gps.latitude * (pi / 180); //assume we have geodetic latitude
double Radlon = gps.longitude * (pi / 180);

double X1= ecef[0];
double Y1= ecef[1];
double Z1= ecef[2];

double c= sqrt(X1*X1+Y1*Y1);
double h= c-Z1/tan(Radlat);

double X2= h*cos(Radlon);
double Y2= h*sin(Radlon);
double Z2= 0.0;

Vec4d ecef_normal(X2-X1,Y2-Y1,Z2-Z1,1.0);

Are my conclusions correct?

What I'm generally trying to do is to determine tilt from gravity normal of some tall object (like Leaning Tower of Pisa) having GPS positions of top point and ground point.


texts:
I am trying to replace the geometry of an autocad mtext object with the insertion point of a block & output it to an ESRI shapefile.  The property fme_geometry for both input mtext & block is fme_point.  BLOCK_INSERTION_POINTS is what the extractor destination and the replacer source is.  Here's the workbench, the mtext is coming from the top, the block from the bottom:

The final shapefile output gives me the correct block insertion points as point features (4 in this case) but without the text attributes I created attached to them, and four more features that have no geometry (fme_geometry = undefined) but do contain the correct text attributes, it looks like this:

The two inspectors in the FME workbench screen capture at top give the following output:

I must be missing something about how GeometryReplacer works or FME is handling the mtext objects, but what?


texts: Road topology problem solution for junction

texts:
I am working with a Road Network. I have digitized the lines with AutoCAD Map 3D. Then Imported it to ArcMap. After entering all kind of database I have created Routing Dataset.  Unfortunately it's not giving the desired route.

The most probable problem may be the junction. I think there are few problems with junction between the road. It's a huge database and it can't be checked manually for every junction. Is there any way to create a network topology to check the junctions easily whether the roads snapped or not?


texts:
What I'd like to know, more than the theoretical limit of the new ArcGIS Pro memory usage, is the real limit.

So, for users of ArcGIS pro, have you run a process using more than 4 GB of RAM, if so, how much? Has anyone experienced any memory limits at all?


texts: Can you synchronize tables between ArcSDE(MSSQL) and PostgreSQL(PostGIS) regularly?

texts:
I have spatial tables in ArcSDE 10.1 (MSSQL 2008 R2) geodatabases and need to update the same spatial tables in PostgreSQL 9.1(PostGIS 2.0.1) geodatabases. Basically I use my ArcSDE environment for editing and geoprocessing,  but my final output needs to be in Postgres(PostGIS) data store for geoserver. The main hurdle for me is keeping my data updated and in sync (no pun intended) between databases.

Using the Data Interoperability extension I set up a spatial ETL model using the "quick export" toolset and can effortlessly move my spatial table(s) from an ArcSDE geodatabase to a PostGIS geodatabase. This is only a one time dump though. How can I  regularly keep my table(s) updated in PostGIS as I make update/insert/delete changes in my ArcSDE geodatabases. I know with ArcGIS 10.1 SP1 you should be able to connect, view and edit spatial data in PostGIS geodatabases if you are using ArcSDE for PostgreSQL. I have ArcSDE 10.1 over MSSQL 2008. Ideally I would need to update my Postgres databases tables once a day to reflect changes made in the ArcSDE geodatabases tables. I am banging my head trying to explore triggers, stored procedures, SQL "magic" and even python scripting.

Please help


texts: Fixing "spatial index grid size is invalid" without Recalculate button

texts:
I've got a python script that is connecting a squence of points to make lines, and putting them into a feature class (based on answers to my previous question). It works wonderfully on a small area that I was using to debug and test (~25 meters, 20 points), but then when I put it into a real-world scenario (~1000 meters, ~500 points), it breaks.

The spatial index grid size is invalid.

Looking into the ArcGIS help suggestions to fix/prevent that, this page suggests recalculating the spatial index grid (directions here), but my polyline feature class doesn't have that option when I am in the "Indexes" tab -- there's no "Recalculate" button at all, and the "Delete" button is grayed out. So how can I fix the error? (I tried simply deleting my original feature class and making a new one, on the theory that my original "little" test case was setting the spatial index grid too small; however, even with a brand-new feature class, I have the same error.)

The problem is more specifically the number of dots I'm trying to connect during one run of the script -- less than ~50 is fine, more than ~100 is not...*

ArcGIS for Desktop with Advanced level license, data stored in feature class in a file geodatabase.


texts: How can I get the middle point between two points

texts:
I am working with a centreline of equidistant points. I want middle points of those points.
For example, if the distance between each of the points is 24.99m. I want to locate the exact middle point at 12.495m. How can I do that with ArcMap tools, not Python?  This is a shapefile.


texts: Metadata (CSW) server + editor alternative to GeoNetwork?

texts:
Can anyone suggest an alternative to GeoNetwork for editing/searching and displaying via a graphical interface?

I quite like it but the "manual update" docs are empty http://geonetwork-opensource.org/manuals/trunk/eng/users/maintainer-guide/updating/manual-update.html and fetching a new 3.2.1 war and pointing it to the old database and copying through trial and error different things from the data folder of the old installation yields blank metadata edit pages, at the best.

Recreating users and metadata permissions and groups and categories isn't so fun so I was wondering if there's an alternative people might recommend?


texts:
I have some cartographic maps I need to fix the position/digitizing.

If you look on the below image you will see that the point layer is not exactly at the center of the map. The exact center is represented by the point shapefile.

Both map and point shapefile are on the same SRC (WGS 84)

The problem is on the map precision.

Is there any way I can easily "drag" the map to center it on the orientation point?

I know I can digitize it again in order to fix the map position, but since I have more than 100 maps with this problem, I am looking for an easier and faster way to it on QGIS?


texts: How do you flip and rotate NetCDF (.nc) files in QGIS?

texts: I am trying to extract Surface PM2.5 values from the NetCDF files available on the Atmospheric Composition Analysis Group's website. However, when I load it into QGIS, it appears to be flipped and rotated. Does anyone know how I can fix this?


texts: I am running a processing PyQGIS script with grass i.pca but I am not getting results in my output folder

texts:
It runs without problem but I don't get the expected pca bands, my output folder is empty, I don't know what to test.
path ='C/Users/Manuel/Documents/INPUT'
output='C/Users/Manuel/Documents/OUTPUT'

bandList = [band for band in
            os.listdir(path) if
            band[-4:]=='.TIF']
input = [os.path.join(path,band)for band in bandList]

def PCA_GRASS(input,output):
    processing.runAndLoadResults("grass7:i.pca",\
    {'-f': False,'n':False,\
    'GRASS_REGION_CELLSIZE_PARAMETER' : 0,\
    'GRASS_REGION_PARAMETER' : None,\
    'input': input,\
    'output': output,\
    'percent': 99,\
    'rescale':[0,400]})
  outList=[band for band in os.listdir(output) if
                      band[-4]=='.tif']
  new_output = [ os.path.join(output,band) for band in
                               outList ]
  for band in new_output:
           iface.addRasterLayer(band,band)

PCA_GRASS(input,output)

When I run the script I open log message and an ogr error appears (can't open _ address of input files_) I don't understand why
(https://i.stack.imgur.com/oAv4s.jpg)
I am working with Landsat 8 band and I have tried other images to see if the code works.  Using processing whit grass I have performed other tasks satisfactorily but this one is causing me many problems.  The images are not the problem when I perform the operation (i.pca) from the processing toolbox, it is done without problem.
I have also changed the syntax with which I enter the input data but nothing.


texts: Distributing aggregated zone data to smaller polygons within these zones, weighted by a field value

texts:
I'm trying to find a hack to better represent the spatial distribution of population and employment density with aggregated zone level data. I'm using ArcGIS Pro (but can use ArcGIS Desktop).

My inputs:

traffic zone (TAZ) polygon shapefile with existing and projected future population/employment
building polygons with heights. I have derived the building volume from this

I would like to distribute the aggregated totals in the traffic zone file to the building polygons, weighted by their volume. I'm struggling to find a way to do this. Anyone have any ideas or resources they can point me to?


texts: Labeling Count values from Arcgis attribute field

texts:
Today at work I was facing some challenges to Label the count values from attribute tables of the feature class. We have all the data stored in geodatabase that is server based. I wrote python code in label expression that did not work out. Challenge is I have a field whose value are repeating,and I want to plot those how many times these values are repeating. Please see below image,

On above image, I want to label A(2), B(2), C(3), D(1). So on the bracket it should show count values.

How to do with the feature class stored in geodatabse that is connection based?


texts: PyQGIS optimize for loop in QGIS script for better performance

texts:
I have a question regarding a for loop that I wrote in a QGIS Python script. The loop works, but it takes too long since I need to let it run over millions of points. So I’m looking for some advice how to improve the coding and hence the performance of the for loop.

The script works as follows: Before the for loop, it takes a point vector layer and writes the points in a point list Pointlist. The for loop over the point list starts: for every point, a viewshed is calculated. The resulting viewshed raster is then loaded in a numpy array and some weighting over clases is done. After that, the sum of the viewshed is calculated. At the end of the loop, the coordinates and the resulting sum are written in a list. The list gets writen outside the for loop (initially the results were written in the result CSV inside the loop, but I thought it is faster outside the loop?).

The script part with the for loop looks like this:

import processing
from qgis.core import *
from qgis.utils import *
from qgis.gui import *
from PyQt5.QtCore import *
import gdalnumeric
import numpy
import csv
import os
import time
import datetime

# viewshed analysis for every point in Pointlist (from Regular Points Grid)
# output is a csv-file with point coordinate and count of viewed pixel

pointgrid_clip = "path_to_my_input_Raster.tif"

inlayer = QgsVectorLayer(pointgrid_clip, "", "ogr")
iter = inlayer.getFeatures()
Pointlist = [f.geometry().asPoint() for f in iter]
totalPoints = len(Pointlist)

Punktliste = []
Coordliste = []
gew_list =[3,2,1,3,2,1]

LoadFile = gdalnumeric.LoadFile
now = datetime.datetime.now

for point in Pointlist:
    n = n-1
    start = now()
    coordstr = '%d,%d' % (point.x(),point.y() )
    extent_viewshed = "%s,%s,%s,%s" % (str(point.x()-(Sichtweite-(pixelWidth/2))), str(point.x()+(Sichtweite-(pixelWidth/2))), str(point.y()-(Sichtweite+(pixelWidth/2))), str(point.y()+(Sichtweite+(pixelWidth/2))))
    viewshed = processing.run("grass7:r.viewshed",
                            {'input':dgm_repclip_ausschluss['OUTPUT'],
                            'coordinates':coordstr,
                            'observer_elevation':Hoehe_Anlage,
                            'target_elevation':1.5,
                            'max_distance':str(Sichtweite),
                            'refraction_coeff':0.14286,
                            'memory':4000,
                            '-c':False,
                            '-r':False,
                            '-b':False,
                            '-e':False,
                            'output':Output_view + "/" + coordstr + "_viewshed.tif",
                            'GRASS_REGION_PARAMETER':extent_viewshed,
                            'GRASS_REGION_CELLSIZE_PARAMETER':0,
                            'GRASS_RASTER_FORMAT_OPT':'',
                            'GRASS_RASTER_FORMAT_META':''})
    raster_file = LoadFile(viewshed['output'])
    werteliste = [w for w in raster_file.flat if not numpy.isnan(w)]
    erg = numpy.histogram(werteliste, bins=[66.1,74.1,82.1,90.1,98.1,106.1,114.1])
    s = (erg[0]*gew_list)
    pixel_count = s.sum()
    Coordliste.append(coordstr)
    Punktliste.append(str(pixel_count))
    ende = now()
    dauer = ende - start
    restdauer = totalPoints * dauer
    feedback.setProgressText('Prozessdauer: ' + str(dauer) + ', verbleibende Zeit~ ' + str(restdauer)  + ', verbleibende Rechnungen: '  + str(n))

os.chdir(Output_Results)
rows = zip(Coordliste,Punktliste)
with open("viewshed_count.csv", 'wt') as csvfile:
    writer = csv.writer(csvfile, delimiter=',', quotechar=',')
    for row in rows:
        writer.writerow(row)


texts: I have a database of raster and vector and I want to use them in my webmapping. There is not a problem to publish the vector using geoserver but I can't publish my raster. Is there possibility to publish postgis raster to geoserver?


texts: How can i set up the "own OSM" for own maps, having no relation to the current real spacial photos

texts:
For some needs we need an OSM-style map displaying the area before the reservoir was constructed and thus displaying the historical shore line. There are some historical aerial photos (not provided by google, but licensed for any usage http://warfly.ru/?lat=55.797473&lon=49.047717&z=12 ) tiled on the manner of google maps, and some old maps stored as large jpeg images (http://xn--80aaaaaa8d3awdlhbbd.xn--p1ai/tatarstan/raritet/volga-kama-pojma-1950).

So, what we should consequently do to be able to produce such kind of online map?

Setup own OSM instance (?) with clear postgis (?) (how?)
Setup the map drawing client that is able to provide tiled WW2 aerial photos and jpeg maps bent to the map? How can we do that, which type of client supports it? If noone, but there are some open source that can be easyly programmed to make this support available with little coding (I hope for my Java coding skills), which opensource client is the most suitable for it?
What is the best way to provide the isohypses information to make the map displaying the  heights? Unfortunately they are not actually stored in osm and will not be available on our sources.
What is the best way to make it available on the www. I know little about postgis, mapnik and tilecache usage, but I'm not sure that it is the best possible technologies stack.

The full guide in answer will be perfect, or even some ideas for some specific points asked are welcome! Thanks in advance!


texts: Geotools always return Multipolygon type for shapefile

texts:
I want to know the true type of Geometry using GeoTools (point, polygon, multipolygon). GeoTools always returns me as a MultiPolygon, and my object has this format for example:

"geometry": {"type": "polygon", "coordinates": [-86.577085000372449, 34.793410999917739], [-86.57692100000412, 34.793837999823893], [-86.576509999936491, 34.794884000396337]

Currently I have the following code:

Map <String, Object> map = new HashMap <> ();
        map.put ("url", file.toURI (). toURL ());

        DataStore dataStore = DataStoreFinder.getDataStore (map);
        String typeName = dataStore.getTypeNames () [0];

        FeatureSource <SimpleFeatureType, SimpleFeature> source = dataStore
                .getFeatureSource (typeName);

        FeatureCollection <SimpleFeatureType, SimpleFeature> collection = source.getFeatures ();
        try (FeatureIterator <SimpleFeature> features = collection.features ()) {
            while (features.hasNext ()) {
                SimpleFeature feature = features.next ();

               feature.getAttribute ("the_geom"). toString ();
}


texts: Calculate Field using Fieldnames(variables?!)

texts:
the Question I have is kind of asked before. But I couldn't figure out the solution.

ArcGIS 10.2.2
Python 2.7.5

I`m adding fields to a table using user input for the fieldnames. Then I  populate those fields and calculate with these fields yet another field.

It looks like this:

fieldname1 = arcpy.GetParameterAsText(2)
fieldname2 = arcpy.GetParameterAsText(3)

arcpy.AddField_management(FC, fieldname1, "SHORT","", "", "",
                          "", "NULLABLE")
arcpy.AddField_management(FC, fieldname2, "SHORT","", "", "",
                          "", "NULLABLE")

In the tool the datatype for the user input is set to "String".

Then I use search and update cursor to populate those fields. Afterwards I want to calculate the next field, so I build a codeblock and an expression and use CalculateField_management:

codeblock = """def myFunction(mthree, mtwo):
    if mthree >= 1 or mtwo >= 3:
        return -3
    elif mthree == 0 and (mtwo==1 or mtwo==2):
    return -2
    elif mthree == 0 and mtwo < 3 and mone >= 5:
    return -2
    if mthree == 0 and mtwo == 0 and mone ==4:
        return -1"""
expression = 'myFunction(fieldname1, fieldname2)'
arcpy.CalculateField_management(FC, "EPK_Gesamt",expression,"PYTHON_9.3",codeblock)

and that doesn't work. I tried several thing such as:

fieldn1 = '"!{}!"'.format(fieldname1)
fieldn2 = '"!{}!"'.format(fieldname2)
expression = 'myFunction(fieldn1, fieldn2)'

or:

expression = "myFunction('!' + fieldname1 + '!','!' + fieldname2 + '!')'

In both cases python returns me text values instead of fieldnames which doesn't work for the calculation in the calculatefield tool.

It works fine if I use the actual fieldnames in the expression variable:

expression = 'myFunction(!example1!, !example2!)'

But of course I want to use the user input, so the user doesnt have to work in the script.
Anyone knows how to deal with that?


texts: Extracting data from attribute table in one shapefile and use it to populate another shapefile attribute table using ArcPy

texts:
I am working on a homework assignment for my coding class, and I am having issues with the last part of the assignment. My job is to create new three shapefiles based on players' positions in a target country. The target country is Sweden, and the positions that we are creating the individual shapefiles from are:

Left-Wing (LW)
Right Wing (RW)
Center (C)

The two fields are the weight and height field, and the weight field is in pounds, and the height field is in inches, and the height field is given in the text in the original shapefile (nhlrosters. hp). The new shapefiles would have the height and weight field, but the height needs to be in cm and kg.
The data for this lesson is located at: https://www.e-education.psu.edu/geog485/node/276
How can I manipulate the string for the height field? Below is the code for the weight field as an example.
import arcpy

fc = r"\\C\\Data":

oldField = "weight"

newField = "weight_kg"

arcpy.AddField_management(fc, newField, "FLOAT")

with arcpy.da.UpdateCursor(fc, (oldField, newField)) as weightCursor:

      for weightRow in weightCursor:

           old_Weight = weightCursor[0]

           new_Weight = old_Weight * 0.453592

           weightRow[1] = new_Weight

           weightCursor.updateRow(weightRow)

For the height column. The height are in feet/inches (example: 6'2")


texts: Is there a way to calculate balances of products between countries?

texts:
I created a dataset with the capacity and production of a product per country.
I am interested in calculating the balances between these different countries.
For example, France has a positive balance of 50, Belgium a balance of 10, and Germany a balance of -50.
I am interested in creating a PyQGIS function that calculates balances between these countries based on the distance to optimize these flows; E.G: Germany to Belgium 10, Germany to France 40, equilibrium reached! Are there any existing functions of which I am not familiar with?

output_folder = 'C:/Users/*/Desktop/Extr-act/Output/'
layer5 = QgsVectorLayer(output_folder + "Distance matrix.shp", "Distance matrix", "ogr")
layer6 = QgsVectorLayer(output_folder + "Centroids.shp", "Centroids", "ogr")

with edit(layer6):
    for f in layer6.getFeatures():
        if f["BC_Balance"] < 0:
            for i in layer5.getFeatures():
                if f["CNTR_CODE"] == i["InputID"]:

                    f["Dest"] = ["InputID"]
                    f["Amount"] = 0
            layer5.updateFeature(f)

I have added the question with the code I came up with so far.

https://drive.google.com/drive/folders/1WCllWC2dgSCsqQ3TokB9obO90jd4ZAoT?usp=sharing


texts:
I couldn't find this question, although I am sure it has been asked.

How could I use a python list to condense this script?  This code block currently works, but is clearly not very efficient.

def Reclass(test):
if test == '2 7':
    return 'Padmounted Single Phase'
elif test == '2 4':
    return 'Padmounted Three Phase'
elif test == '2 2':
    return 'Padmounted Three Phase'
elif test == '2 1':
    return 'Padmounted Three Phase'
elif test == '1 4':
    return 'OH Single Phase'
elif test == '1 2':
    return 'OH Single Phase'
elif test == '1 1':
    return 'OH Single Phase'
elif test == '1 5':
    return 'OH Second Phase'
elif test == '1 3':
    return 'OH Second Phase'
elif test == '1 7':
    return 'OH Third Phase'

I was attempting to use (as an example of a list):

elif test == ['1 1', '1 2', '1 3']

but the code was being skipped over and not working.


texts:
I have downloaded MODIS data (from ladsweb.modaps.eosdis.nasa.gov) which I would like to view in ArcMap 10.4.  It's MOD14 thermal anomalies data.

Currently I can place the data on a map, but it appears to be in the middle of the Pacific Ocean (which is not the correct placement of my data). This suggested to me that the coordinate system or projection was incorrect. However, I'm not able to determine what the geographic coordinate system should be.

The spatial extent should cover Kansas:

I've tried several things:

I went back to the NASA page where I downloaded the data. It does not document the coordinate system being used.
I searched for coordinate systems containing the word "MODIS" in ArcMap. I also used the search tool to look for other MODIS-related subjects in Esri's help.
I googled my problem. Many pages suggest using R scripts or third-party tools to transform the data prior to being used in ArcMap. I'm not sure what this is about or why it's necessary - it seems like a lot just to import some data.
I recently went through several of ESRI's online training modules about using geographic coordinate systems and projects, but didn't find anything helpful.

Do I really need a third-party tool or script to use MODIS data in ArcMap? If not, what do I need to do to work with MODIS in ArcMap?


texts:
I understand this question has been asked a few times but I am having trouble with this issue in QGIS 3.0.3.

These are the current CRS of my layers:

Open Street Map (XYZ Tile) - EPSG: 3857 - WGS 84/Pseudo Mercator
Addresses - EPSG: 4326 - WGS 84
Road Data - EPSG: 4326 - WGS 84

I am trying to calculate distances using NNJoin tool from addresses to the road in the state of Missouri. I need my distance in meters but it only shows up in decimal degrees. I have enabled "on the fly" projection, changing CRS to a UTM zone (central Missouri), and making sure all the settings are presenting distances in meters (see image attached). Yet, I continuously get the warning that the measurement for the NNJoin tool will be in decimal degrees.


texts: What spatially-aware GUI administration softwares exist for PostgreSQL/PostGIS?

texts:
Most of my experience with PostgreSQL has been pgAdminIII and psql, but I'm interested to know what others use.

Is there PostgreSQL equivalent to SQL Server Management Studio (SSMS) and its ability to display, in a limited way, results containing geometry (free or paid)?

Or do most people just run a GIS client on the side for this purpose, as this question seems to indicate: GUI tools for PostgreSQL/PostGIS??

I looked down through the PostgreSQL software catalog, but I didn't find anything (other than broken links).


texts: How to interpret the Maximum flow path length raster, from the SAGA Maximum flow path length tool

texts:
Does anyone know how to interpret the Maximum flow path length raster, after running the SAGA Maximum flow path length tool?
I have used this tool and received a raster but I'm not sure how I can get the maximum flow length (possibly as a line or even just a number) for this catchment. The raster values vary from 0 to 139. But I'm sure the flow path should be of order of 600-700 feet. I attached the output result if it helps.
As input, I just had the filled elevation raster.


texts: How to automate writer output directory and filename in FME Desktop?

texts:
I built a workspace in FME Desktop 2017 for re-projecting DWG files between different coordinate systems and I need to automate the writer output folder and filename of the dwg in FME desktop.

I need to compose automatically the writer output path using: 'reader input path/myfolder/reader filename.dwg'

Is there anyway to do this in FME Desktop?


texts: CartoType - set min/max on map view while Panning

texts:
Working with the .NET CartoType SDK on Windows. How can I restrict panning? In the "Windows Maps" app I am able to pan outside of the map boundaries. I would like to prevent this.

The C# demo program supplied with CartoType contains this code to pan the map in response to a mouse drag event, and I've added the code starting with the comment "clamp view to extent of map", but the map is still panned too far. The aim is to prevent any edge of the display going beyond the edge of the map.

private void Form1_MouseMove(object sender, MouseEventArgs e)
    {
        if (e.Button == MouseButtons.Left)
        {
            m_map_drag_offset_x = e.X - m_map_drag_anchor_x;
            m_map_drag_offset_y = e.Y - m_map_drag_anchor_y;

            if (m_graphics_acceleration)
            {
                m_framework.Pan(-m_map_drag_offset_x, -m_map_drag_offset_y);
                m_map_drag_offset_x = 0;
                m_map_drag_offset_y = 0;
                m_map_drag_anchor_x = e.X;
                m_map_drag_anchor_y = e.Y;

                // clamp view to extent of map
                var state = m_framework.GetViewState();
                var map_extent = new CartoType.Rect();
                m_framework.GetMapExtent(map_extent, CartoType.CoordType.Degree);
                var x = state.m_view_center_degrees_x;
                var y = state.m_view_center_degrees_y;
                if (x < map_extent.MinX)
                    x = map_extent.MinX;
                else if (x > map_extent.MaxX)
                    x = map_extent.MaxX;
                if (y < map_extent.MinY)
                    y = map_extent.MinY;
                else if (y > map_extent.MaxY)
                    y = map_extent.MaxY;
                m_framework.SetViewCenter(x, y, CartoType.CoordType.Degree);

            }
            else
                Draw(m_map_drag_graphics);
        }
    }


texts: Using R to plot maps using different projections

texts:
I currently plot global maps in R using the standard projection in maps (). As far as I understand this is the mercator projection.
Now I want to use a projection that better represents the areas of different countries. I tried:

map('world',proj="mollweide")

But then I cannot plot the shapefile I used to overlay with the initial map or the points from my data.frame.
Do I need to change the projection of everything I want to plot on top of the global map which does not uses the default projection? If so, how do I do that?
Is maps() the best function for this purpose?


texts: Differentiate sectors with function editor of QGIS

texts:
I am new to QGIS.

I have a layer that, among other fields, has two with the structure below:

id | Class
1  | 2
2  | 1
3  | 2
4  | 5
5  | 6
6  | 2
7  | 3
8  | 4

The id column represents specific regions of a layer and class has values ranging from 1 to 7. I would like to differentiate the regions based on the value of the layer column. I am trying to do it through the function editor, but I am not able to reference the columns in the code.

How do I do this task?


texts: Connect points to a line with shortest distance using QGIS?

texts:
On the picture below, I want to let my point connect to the line with the shortest distance and show the distance to me.

I tried NNjoin, MMQGIS, v.distance But I still can't get the distance.

How do I connect them and get the shortest distance?

After that, I used the field calculator, but I think the distance is wrong

I measured it by myself, is any details can tell me to calculate the distance?


texts: PostGIS: intersection of a polygon with 2 rasters

texts:
I'm trying to get the intersection of a Polygon with 2 raster layers in PostGIS.

The rasters have the same size and pixel size and I'd like to retrieve for the polygon a list of something like that:

1; 2
1; 2
1; 3
2; 2
...

Where each line is a pixel intersected by the polygon in both layers and the values is a combination of the pixel value from raster 1 and the pixel value from raster 2.

I don't really know what kind of request I have to do to achieve that...

EDIT: Maybe this query can answer my question:

SELECT
    ST_DumpValues(ST_Clip(nh.rast, 1, geom, true), 1, false),
    ST_DumpValues(ST_Clip(nt.rast, 1, geom, true), 1, false)
FROM
    niv_hydrique AS nh,
    niv_trophique AS nt,
    ST_Transform(ST_GeomFromText('%s', %s), 31370) AS geom
WHERE
    ST_Intersects(nh.rast, geom) AND ST_Intersects(nt.rast, geom)

Where niv_hydrique, niv_trophique are my 2 rasters and I construct my polygon from a WKT representation.

I can then programmatically extract the values for each pixel in the two 2-dim arrays returned.


texts: How to export georeferenced Raster layer using QGIS2Leaf Plugin?

texts:
I got excellent results while exporting vector layer with styling through qgis2leaf plugin. Whereas when I try to do the same with a geo-referenced raster file (geotiff) which has a projection of EPSG:4326, I am getting a Python Error. I tried reprojecting it to EPSG:3857, even tried changing the format of raster file, but getting the same error. here is a screenshot of the error

Is there any additional GDAL Framework which needs to be installed / updated? Any guidance would be of great help as i am trying to have a transparent existing land use map over OSM to check disparities.

Abhijit


texts: Supporting NON GIS capable people to draw/make maps?

texts:
As GIS/cartography capable people... what do you advise when faced with "can you show me how to make a simple map?" questions from non-GIS capable people? Some answers I can think of (depending on the enquirer) are:

Go to Google Maps and use their tools, then print the page (but take note of copyright issues).
Export a raster image from Openstreetmap (using the Bigmap service or the openstreetmap.org share dialogue), display this in Powerpoint as an image and draw on top of it using Powerpoint tools.
Do this but put the image into your favourite powerful raster graphics application to edit.
Export a vector image from Openstreetmap (using the openstreetmap.org share dialogue) and import this into your favourite vector graphics app so you can draw on it.

I'm looking for answers which do NOT involve downloading new software or setting up a simple GIS  system or learning a new complicated process. I can think of solutions using web-GIS tools or QGIS or specially designed systems - but in this case I'm looking for answers which would apply to any random enquirer who has basic computer skills and access to relatively standard software with which they are already familiar.

It might not be entirely obvious but I'd also like comments about simple common pieces of software which work well for drawing markers/lines/words/etc on top of either a raster or vector map image - or an equivalent process with which a map can be customised. I mention this because I recently discovered how useful Powerpoint and its equivalents can be (this is not obvious) for handling documents which are heavily based on annotated photographs (but which are never intended for use as a presentation).


texts: Feature space selection for river sand classification using multispectral imagery

texts: I want to perform supervised classification for river sand detection. I have 3 m spatial resolution Surface Reflectance multispectral imagery with B, G, R and NIR bands (Planet data). Which all features should I choose for classification input that would give best separation for my sample dataset?


texts: Setting ArcMap to get control points in pixels?

texts:
I need to export control points from ArcMap to use it with Matlab. I want to get these points in pixel not in the meter as you can see the below image

these images are 1400X900 px. I expect point 1 = (0,0) and 3 (900,1400).

what should I do to make the control points in pixel?


texts: T-rex tile server with PostGIS (Windows) : utf-8 encoding error

texts:
I want to use T-rex as vector tileserver on Windows.
I load osm.pbf file (from geofabrik) to PostGIS database with osm2pgsql successfully.
When I run the server with the following commands, I get the error message :

Command :

t_rex serve --dbconn postgresql://Myuser:Mypass@localhost/Mydb

Error :

ERROR IO error: invalid utf-8 sequence of 1 bytes from index 34
Even when I just want to generate the config file, same error :
t_rex genconfig --dbconn postgresql://Myuser:Mypass@localhost/Mydb
I suspect this is a postgis setup problem ... that may not have anything to do with T-rex?
Unless it is an import configuration problem?
Note that using T-rex directly with a pbf file does not pose this problem (t_rex serve --datasource My.osm.pbf).


texts:
I was just curious as to whats the best way within a python script to get the output a dbf file/ table while the input is a shapefile. (if possible which i'm assuming it is)

Is there a way of converting a shape file into a table, or by selecting the attribute table of the shapefile and  just outputting that.

Or am i going to need to create a table and then import the values?


texts: ArcPy select by location to copy attribute values from one shapefile to another

texts:
I have two shapefiles, lines and points.  I want to copy over the lines number to the points that are near each line.

I can use select by location in ArcMap and label manually for a small scale, but need to automate on a larger scale.


texts: Transfer fields from one layer to another ArcGIS Pro SDK

texts:
I am creating a new layer and I would like for the new layer to take on the fields schema from the input layer. How can I accomplish this in ArcGIS Pro SDK C#?

This is the code that I have to retrieve the fields for the inRoutes. How do I apply these fields to the newly created feature class?

var fields = inRoutes.GetFeatureClass().GetDefinition().GetFields();
var newfc = await Geoprocessing.ExecuteToolAsync("CreateFeatureclass_management", Geoprocessing.MakeValueArray(fcPath, filename, "POLYLINE", "", "ENABLED", "DISABLED", sr), null, null, null, GPExecuteToolFlags.AddOutputsToMap);


texts: How can I calculate the average speed in various legs of a bus route?

texts:
I'm working on analyzing NYC MTA Bus traffic data.

There are numerous buses traveling along the route, emitting gps data at random locations.

I have about 30k records with bus_line, vehicle_id, longitude, latitude, bus_speed and timestamp. I refer to this dataset as dataset1

Calculating the speeds for an individual bus is trivial but I can't think of a simple way of averaging the speed in a particular leg along the route.

One idea I'm thinking about is to create another dataset (i.e., "dataset2") with longitude & latitude in 50 ft intervals. Then find a way to assign its points to the points from dataset1 based on proximity. At this point, it would be possible to find the average speed for each point because there will be multiple records with the same longitude and latitude.

The biggest challenge I see is in assigning the points from dataset2 to the points in dataset1: you need to sort the points exactly in the order they "appear" in the route. You can't simply order by longitude or latitude. You can sort conditionally such that if the bus is going west, you sort by longitude. If it's going north, you sort by latitude. However, I'd like something simpler if it exists.
I'd appreciate any help!

Thanks.


texts: error "ValueError: API 'QDate' has already been set to version 1" when converting to standalone qgis script

texts:
I'm trying to convert a python script for qgis (2.18.2 Las Palmas) to executable (.exe) using pyinstaller

This is a very simple example (I'm just a beginner):

from qgis.core import *
import os
print "prova"

this is the coomand I use to convert esempio.py to esempio.exe:

pyinstaller esempio.py --onefile

when I run esempio.exe the following message appears:

Traceback (most recent call last):
  File "esempio.py", line 1, in <module>
  File "c:\progra~1\qgis2~1.18\apps\python27\lib\site-packages\PyInstaller-3.4-py2.7.egg\PyInstaller\loader\pyimod03_importers.py", line 395, in load_module
    exec(bytecode, module.__dict__)
  File "qgis\__init__.py", line 26, in <module>
  File "c:\progra~1\qgis2~1.18\apps\python27\lib\site-packages\PyInstaller-3.4-py2.7.egg\PyInstaller\loader\pyimod03_importers.py", line 395, in load_module
    exec(bytecode, module.__dict__)
  File "qgis\PyQt\QtCore.py", line 28, in <module>
ValueError: API 'QDate' has already been set to version 1
[10208] Failed to execute script esempio

Can someone help me


texts: Removing Esri logo from ArcGIS Viewer for Flex 2.2?

texts:
How can we remove the Esri logo from ArcGIS Viewer for Flex 2.2?

Logo layer is set to be visible by default.

Is there a way to turn it off?


texts: How to sort attribute table using ArcObjects/VBA code?

texts:
I'm using ArcMap 9.3.

Shift the rows order display (ascending/descending) for a given field is as simple as double clicking on the screen.

However I'd like to know how to perform this simple operation through automation, in my case, using VBA code.

I've done searches over several interfaces such as ITableDefinition, ITableSort, IIndex with  no success.  ITableSort does exactly this but I did not find out an interface to pass its arguments to IDisplayTable.

Can someone help me find the parameters to pass?


texts: Get latest tiles from TileMill's tile server in Leaflet

texts:
My leaflet app gets the tiles from TileMill's tile server in the standard way:

L.tileLayer('http://<myip>/tile/SteveCountryVic/{z}/{x}/{y}.png', { ... });

These tiles are often old, not showing the latest styles. I'm not sure if they're coming from TileMill itself or some intermediary cache.

I notice that TileMill includes an 'updated' parameter:

http://<myip>/tile/SteveCountryVic/{z}/{x}/{y}.png?updated=1372394902000

Manually copying that updated timestamp means I get the latest timestamp. But what's a good way to not have to manually update the URL every time I change the style?


texts: Rotate marker symbols individually in QGIS

texts:
Is it possible to rotate marker symbols (e.g. a triangle) in QGIS individually with an attribute field like in ArcGIS?

Edit: I would like to rotate each marker symbol individually in category A with an attribute field X, where the angle is defined.

This is how it is done in ArcGIS:
"To rotate marker symbols individually based on the value of a field in the point layer's attribute table click Layer Properties dialog box > Symbology tab > Advanced > Rotation. On the Rotate dialog box, the Rotate Points by Angle in this field specifies the field that holds angular values and the Rotation Style controls the origin and direction of rotation" (ESRI Help).


texts: Has "add polygon attributes to points" been removed from QGIS?

texts: I've been using "add polygon attributes to points" in qgis for a while, but hadn't needed it for a while (@Underdark provided this); in 2.13.0 master on ubuntu it seems to be missing from the processing toolbox. Am I missing something? Do I need to install something special?


texts:
I have an image where pixel values are integers, they go from 0 to 6 and I also have some 999. I want the value 999 to be transparent, which I could achieve by myself. I would like values 0 to 6 to be different colors. I have no idea how to do this.

Can anyone help me?


texts: S-57 import in QGIS; selecting one layer from multiple files?

texts:
I am building a nautical basemap that is build out of multiple S-57 cells (*.000 files).

I go through to following steps:

I download the files from a server and copy all the files (65 files) in a folder
I open QGIS and go to folder in "Panel Browser"
I select all the files and drag them to my "Panal Layers"
I pop-up opens "Selecteer toe te voegen vectorlagen" (that should translate to: Select vectorlayers to be opened"
I manually, handselect, the layer i want, which is: DEPARE.
After the 1st popup, 64 more automaticly open, so I repeat the process 64 times.
After all this is done. I have all my DEPARE layers on screen.
I select all the DEPARE layers (vectorlayers) and merge them to one shapefile.

Is there a way in QGIS to prevent step 6 of the above?

I would like to add all the files to my project but have QGIS automaticly pull the DEPARE layer from all the files, so I don't have to handselect them 65 times.


texts: Change layer styles based on ComboBox search?

texts:
I'm using a combination of OpenGeo Server components and hand written GeoExt code to create a map interface. I would like to have a ComboBox that queries the OpenGeo PostGIS database in order to create a list of counties in a state. As the user begins typing the ComboBox searches the store of counties to narrow down to the one they are looking for. Upon selecting the county I would like to change the style of the layer in order to show the selected county.

Using the GeoExplorer I am able to create various styles in a single layer that relate to certain rules. For example:

I have a layer that references 100 counties in a state. The layer is made up of 100 different shape geometries, each of which defines a single county. I have one layer style that outlines the counties like this:

<sld:NamedLayer>
<sld:Name>opengeo_CountiesNC_OSM</sld:Name>
<sld:UserStyle>
  <sld:Name>opengeo_CountiesNC_OSM</sld:Name>
  <sld:Title>A transparent stroked polygon style</sld:Title>
  <sld:IsDefault>1</sld:IsDefault>

I have another layer style that refers to only a specific county and outlines it. This style uses a rule "Match: Name = Durham" in order to determine which county to outline.

<sld:NamedLayer>
<sld:Name>County Lines_5f024dcb</sld:Name>
<sld:UserStyle>
  <sld:Name>County Lines_5f024dcb</sld:Name>
  <sld:Title>Durham</sld:Title>
  <sld:Abstract>County lines for only Durham County</sld:Abstract>
  <sld:FeatureTypeStyle>
    <sld:Name>name</sld:Name>

I would like to have a style that represents each of the counties in the layer so that they could be individually selected and shown on the map. Rather than have 100 different layers I believe it is easier to have 100 styles.

Is there some way to change which style of a layer is being displayed via GeoExt or OpenGeo JavaScript code? Having a listener on the ComboBox to read the selected record and display the corresponding layer seems ideal. I'm not sure what access I have to the layer styles that are kept in the GeoServer from the web application side of things. What calls can I make in JS to appropriately retrieve and display data?

Edit
When creating a WMS Layer there is an optional param "styles" that allows for an array of String style names:

new OpenLayers.Layer.WMS("County Lines",
                "geoserver/wms",
                {layers: "opengeo:County Lines",
                    transparent: true,
                    styles:[
                        "opengeo_CountiesNC_OSM",
                        "County Lines_Durham"
                    ]},
                {
                    visibility: true,
                    displayInLayerSwitcher: false,
                    opacity: .7
                }
            )

With multiple styles imported it should be possible to switch between them?

EDIT 2

Ended up doing something similar to the WMS Store. Because I already have a Google Geocoder I discovered that I could hook into the results of that ComboBox to determine when a county had been search for. Once a search had been performed the results can be trimmed and inserted into a merge request.

if(str.match(/(\D\s)+County, NC, USA/))
                    {
                        console.log("Match! It's a county");
                        console.log("County Lines_" + str.substring(0, str.indexOf('C')));
                        mapPanel.map.layers[2].mergeNewParams({styles: "County Lines_" + str.substring(0, str.indexOf('C'))});}


texts:
While trying to read a layer from PostgreSQL, in order to select one polygon and zoom to this feature, I see the entire layer but not the selected feature.
I am using an unconventional way solve the problem - I am running the three lines I marked as comments in my example.

Can anyone offer me a better way?

I am working with Windows 10, QGIS 2.18.22, PostgreSQL 10.

This is my code that I am running in Python Consol Editor:

from qgis.core import *
from console import console

uri = QgsDataSourceURI()
uri.setConnection ("127.0.0.1", "5432", "DB_Layers", "postgres", "myPassword")

uri.setDataSource("public",  'blocks_lay', "wkb_geometry", "")

layer1 = QgsVectorLayer ( uri.uri(False), "blocks", "postgres")
QgsMapLayerRegistry.instance().addMapLayer(layer1, True)

crs = QgsCoordinateReferenceSystem(2039)
iface.mapCanvas().mapSettings().setDestinationCrs(crs)

layer1.setCrs(crs)
expressionCur = 'block_num = 6018'

it = layer1.getFeatures(QgsFeatureRequest().setFilterExpression(expressionCur))
layer1.setSelectedFeatures ( [ f.id() for f in it ])

boxWin = layer1.boundingBoxOfSelected()
# console.show_console()
# print 'aaaa'
# console.show_console()
iface.mapCanvas().setExtent(boxWin)

iface.mapCanvas().refresh()


texts:
Context

ArcGIS 9.2 Build 1500

ArcObjects Windows Forms application
running outside of ArcMap, similar to ArcEngine product.
Accessing mxd to print out page
layout to A3 paper.  Mxd is saved
with same settings as printer, see
below:
Original code was copied from PrintActive View sample on ESRI site.  This sample does run in the ArcMap environment, but I don't think this should have a bearing on my problem.

Actual Problem

For some strange reason when my code prints out the map - it always prints to B5 which is obviously wrong as it should be A3.  Printer complains that it doesn't have B5 - no map gets printed and misery ensues!
On further inspection of the code it seems that something strange is happening, the code below is a cut down version of the sample code just to show how you get the problem:

IPrinter docPrinter;
IPaper docPaper;

/* printdocument is from the .NET assembly system.drawing.printing */
System.Drawing.Printing.PrintDocument sysPrintDocumentDocument;

/* Now we need to get the default printer name.  Since this is a generic command,
 * we can't use the printername property of the document.  So instead, we use the
 * System.Drawing.Printing objects to find the default printer.
 */
docPrinter = new EmfPrinterClass();
sysPrintDocumentDocument = new System.Drawing.Printing.PrintDocument();
docPaper = new PaperClass();

/* testing to see if printer instantiated in sysPrintDocumentDocument is the
 * default printer.  It SHOULD be, but this is just a reality check.
 */
bool isDefault = sysPrintDocumentDocument.PrinterSettings.IsDefaultPrinter;

if (isDefault)
{
    //FIRST STEP
    //Set docPaper's printername to the printername of the default printer
    docPaper.PrinterName = sysPrintDocumentDocument.PrinterSettings.PrinterName;

    //Need to have the same paper form as the mxd
    docPaper.FormID = (short) mapDocument.PageLayout.Page.FormID;
}

/* Now assign docPrinter the paper and with it the printername.  This process is two steps
 * because you cannot change an IPrinter's printer except by passing it as a part of
 * the IPaper.  That's why we setup docPrinter.Paper.PrinterName first.
*/
//SECOND STEP
docPrinter.Paper = docPaper;

when the code sets the docPaper.FormID (it's defaulted to 1 which is esrPageFormLetter) and overwrites this value with the value in the mxd (IMapDocument) mapDocument.PageLayout.Page.FormID - in this case its 13 (esrPageFormSameAsPrinter), that makes sense as we saved the mxd to "Use Printer Paper Settings" (see image above).

The 13 value that is passed seems to cause some conflict in translation to the IPaper object - because it was passed 13 it thinks that is B5 and not to try and infer out the printer settings which should be A3.

The best I can make out of this - is that it must have something to do with Windows paper size constants e.g. DMPAPER_B5 - 13 -  B5 (JIS) 182 x 257 mm.

Does anyone have a clue how to get round this problem?


texts: Using onRectangle for user input in ArcMap/ESRI addin

texts:
In an effort to try to add user input to my project, I'm developing an esri python addin to generate rectangles to my map.  I'm on step one now, just exploring the addin functionality.

I tried making an addin 'tool' precisely per the ESRI sample code, but the button isn't working.  It appears on my ArcMap desktop, but when I click it, nothing happens.

I just really want to know what's happening with the onRectangle, and how I can get it to accept user input.

Here's the sample code that I essentially stole from ESRI:

import arcpy
import pythonaddins

class fishnetsExample(object):
    def __init__(self):
        self.enabled = True
        self.cursor = 3
        self.shape = 'Rectangle'

    def onRectangle(self, rectangle_geometry):
        """Occurs when the rectangle is drawn and the mouse button is released.
        The rectangle is a extent object."""
        extent = rectangle_geometry
        # Create a fishnet with 10 rows and 10 columns.
        if arcpy.Exists(r'in_memory\fishnet'):
            arcpy.Delete_management(r'in_memory\fishnet')
        fishnet = arcpy.CreateFishnet_management(r'in_memory\fishnet',
                            '%f %f' %(extent.XMin, extent.YMin),
                            '%f %f' %(extent.XMin, extent.YMax),
                            0, 0, 10, 10,
                            '%f %f' %(extent.XMax, extent.YMax),'NO_LABELS',
                            '%f %f %f %f' %(extent.XMin, extent.YMin, extent.XMax, extent.YMax), 'POLYGON')
        arcpy.RefreshActiveView()
        return fishnet


texts: Classifying a GeoJSON LineString type file based on Strahler Order with Python

texts:
I am working on a Python project working with a GeoJSON file (Linestring type).
The aim of this project is to classify lines based on Strahler Order (Wikipedia : https://en.wikipedia.org/wiki/Strahler_number). Here below a picture to understand it.

I run a Python script that I wrote to detect OERDER 1 but when I test it on a line that I know it got Order 1, I always get 0. I don't know where is the problem in my code.
Here below my code :
import json

file = open('C:/Users/LENOVO/Desktop/data_with_subgraphs5.json')
data = json.load(file)

def search(data):

    featureA = []
    featureB = []

    for feature in data['features']:

    featureA = feature['geometry']['coordinates'][0]

        for ftr in data['features']:

           featureB = ftr['geometry']['coordinates'][0]

           if featureA[0] in featureB or featureA[-1] in featureB :

               feature['properties']['Strahler'] = 1

           else :

               feature['properties']['Strahler'] = 0

    return data['features'][372]['properties']


texts:
the labels in GeoNetwork are defined as:

<element name="gmd:codeSpace" context="gmd:RD_Identifier">
  <label>Namespace</label>
</element>

but the problem is we also have a Regulation Source which is saved in the metadata in a gmd:codeSpace tag and if we create a new element for that it gets ignored since both have the same name and the same context (gmd:RD_Identifier).

Is there a way to work around this?


texts: Does anyone have any tips on figuring out why a Vector Layer didn't load in Pyqgis?

texts:
This is a general question.  I'm stepping through the pyqgis code and the layer isValid() is always False.  I don't see any error messages showing any reason why it wouldn't load.  Other similars have loaded without problems.

Just trying to see if I need to actually download the QGIS source and start stepping through that as well...?  (i hope not)

the first layer loads without problems.  it's referencing a postgis table.  the second doesn't but the another_table is a view. perhaps that makes a difference but there are no error messages or exceptions letting me know what's wrong...

Nothing fancy:

    #first vLayer
    self.uri.setDataSource('public','some_table','location','')
    self.first_layer = QgsVectorLayer(self.uri.uri(), 'Some Locations', "postgres")

    QgsMapLayerRegistry.instance().addMapLayer(
        self.first_layer)

    #second Layer
    self.uri.setDataSource('public', 'another_table', 'location','')
    self.second_layer = QgsVectorLayer(self.uri.uri(), 'Another set of Locations', "postgres")

    QgsMapLayerRegistry.instance().addMapLayer(
        self.second_layer)


texts: Why does Google Maps display strange, jagged state boundaries in the Great Lakes?

texts:
When I search for a state on Google Maps, I am shown the state outlined in red. I've noticed that the states that border the Great Lakes have boundaries that protrude into the lakes in strange ways.
Here are some examples:

What is going on here?


texts:
My organisation uses a SAAS and part of that service publishes data via a WFS using Geoserver 2.7.0. I have sent a GetFeature request via a browser (IE11) and it returns the GML however, if I refresh the link the FID changes:

:wfs_BBQ gml:id="wfs_BBQ.fid-74424e53_159f7837ba8_350">

:wfs_BBQ gml:id="wfs_BBQ.fid-74424e53_159f7837ba8_ec3">

I am using QGIS 2.18.3 (running on Win 10) and I'm told QGIS needs persistent FIDs (see: Duplicate WFS Featues in QGIS 2.18.3) and that Geoserver should create a persistent FID if the WFS has a Primary Key. The data does have a Unique ID (integer) so I'm hoping this meets the requirement.

How do I ensure the FIDs are persistent?


texts: Classification: Layer error: Classifier training failed: 'Only one class.'

texts:
i am new in GEE and try a tutorial from this https://www.gears-lab.com/intro_rs_lab4/ but the result show error classifier training failed: 'Only one class.'. So i try to fix it and combine with this tutorial https://developers.google.com/earth-engine/classification but again the result show error training failed: 'Only one class
this is my code https://code.earthengine.google.com/5412d64d8ba217f879be98f6ba491146
im stuck.
sorry my bad english.


texts:
I usually create hexbins using the MMQGIS tools in QGIS. I would create a hex layer and intersect it with a shapefile. I've tried to create a hex shapefile and importing it into a PostGIS and tried to do an ST_Intersects but the result usually come up empty.

I've tried looking at other tutorials on hexbinning on PostGIS like this and this but they create a blanket hex bins that isn't .

Is there a way to create hex bin from PostGIS over a shapefile (truncation of hexes are ok).


texts: How to make a match the GPS location to a model?

texts:
So I have a map like this where I need to match the GPS position of a device to a specific position in a model:

This is roughly the model of the place where I want to use the GPS, is there some way I can achieve what I want?
BTW this is the place where I want to use the GPS of the device

I'm not sure if this is the right place to ask this question but if not, please point me to it.


texts: QTableView / QTableWidget alternative for floats

texts:
Writing a QGIS plugin with PyQGIS I need a table to store and view values that are floats. I found only these two models (QTableView / QTableWidget), but they only store QString type of data, and therefore sorting of numbers gets wrong (eg. 10, 23512, 24, 35, 6874 would be correct order while it's obviously not).
Main problems that I have are:

Proper sorting is very important
Size of values stored can be any (eg. can't say it will be always 4 digits with 2 precision digits)

Is there any alternative to store numbers in tables in QGIS plugins?


texts: How to import and layer stack seperate ndvi bands to grass gis?

texts: I have 20 separate single band NDVI files. I have to import this into GRASS and stack them  into a single 20 band image.  How can I do that?