texts: Mixed geometry types styling with GeoServer 2.2

texts:
I have a PostGIS table with geometry column. In this column I have geometries of several types (polygons, points, lines). I publish this this table on GeoServer.

Now I want to create OpenLayers layer with this data.  But how to make a style for this data?

I've already read http://docs.geoserver.org/stable/en/user/styling/sld-tipstricks/mixed-geometries.html  and don't want to create gtype in PostGIS if this possible.

I try to create new style within GeoServer but I think that do something wrong:

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
 xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
 xmlns="http://www.opengis.net/sld"
 xmlns:ogc="http://www.opengis.net/ogc"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
<!-- a Named Layer is the basic building block of an SLD document -->
    <Rule>
     <ogc:PropertyIsEqualTo>
      <ogc:Function name="geometryType">
       <ogc:PropertyName>the_geom</ogc:PropertyName>
      </ogc:Function>
      <ogc:Literal>Polygon</ogc:Literal>
     </ogc:PropertyIsEqualTo>
     <PolygonSymbolizer>
       <Stroke>
          <CssParameter name="stroke">#0000ff</CssParameter>
          <CssParameter name="stroke-width">1</CssParameter>
        </Stroke>
     </PolygonSymbolizer>
    </Rule>
</StyledLayerDescriptor>

And get a error:

line 10: cvc-complex-type.2.4.a: Invalid content was found starting with element 'Rule'. One of '{"http://www.opengis.net/sld":Name, "http://www.opengis.net/sld":Title, "http://www.opengis.net/sld":Abstract, "http://www.opengis.net/sld":NamedLayer, "http://www.opengis.net/sld":UserLayer}' is expected.
line 11: cvc-complex-type.2.4.a: Invalid content was found starting with element 'ogc:PropertyIsEqualTo'. One of '{"http://www.opengis.net/sld":Name, "http://www.opengis.net/sld":Title, "http://www.opengis.net/sld":Abstract, "http://www.opengis.net/sld":LegendGraphic, "http://www.opengis.net/ogc":Filter, "http://www.opengis.net/sld":ElseFilter, "http://www.opengis.net/sld":MinScaleDenominator, "http://www.opengis.net/sld":MaxScaleDenominator, "http://www.opengis.net/sld":Symbolizer}' is expected.

What am I doing wrong?

UPDATE

I try to show a points and make this rule for this geom type:

 <Rule>
        <ogc:Filter>
            <ogc:PropertyIsEqualTo>
                <ogc:Function name="geometryType">
                    <ogc:PropertyName>the_geom</ogc:PropertyName>
                </ogc:Function>
                <ogc:Literal>Point</ogc:Literal>
            </ogc:PropertyIsEqualTo>
        </ogc:Filter>
        <PolygonSymbolizer>
        <Stroke>
            <CssParameter name="stroke">#0000ff</CssParameter>
            <CssParameter name="stroke-width">2</CssParameter>
        </Stroke>
        </PolygonSymbolizer>
    </Rule>

But not see then on my map. What's wrong?


texts: Joining 2 tables to create layer in Geoserver with MultiPolygon and Point data types

texts:
I got problem when using the SQL View function. My objective is to create layer in geoserver from joining 2 tables.

I have 2 tables, table_a:

                                 Table "public.table_a"
         Column          |             Type             |        Modifiers
-------------------------+------------------------------+------------------------
 id                      | integer                      | not null default
 name                    | character varying(30)        | not null
 remarks                 | text                         |
 boundary                | geography(MultiPolygon,4326) |

And table_b:

                           Table "public.table_b"
         Column          |         Type          |           Modifiers
-------------------------+-----------------------+-------------------------------
 id                      | integer               | not null default
 sequence                | integer               | not null
 point                   | geography(Point,4326) |
 table_a_id              | integer               | not null

Here is the query for the SQL View:

SELECT
    table_a.name, table_a.remarks, table_a.boundary,
    table_b.sequence, table_b.point, table_b.table_a_id
FROM
    table_a
        INNER JOIN table_b ON (table_a.id = table_b.table_a_id)

And the result:

      Name         | Type
-------------------+-----------
table_a.name       | String
table_a.remarks    | String
table_b.sequence   | Integer
table_b.table_a_id | Integer

The problem with the result is the table_a.boundary column and table_b.point column cannot be display.

The result is different when I tried with Publish function for each table:

Publish table_a:

      Name         | Type
-------------------+-----------
name               | String
remarks            | String
boundary           | MultiPolygon

Publish table_b:

  Name             | Type
-------------------+-----------
sequence           | Integer
Point              | Point
table_a_id         | Integer

My Question is:

What is the query to create the MultiPolygon data type in GeoServer SQL View from geography(MultiPolygon) in postgres data type?
What is the query to create the Point data type in GeoServer SQL View from geography(Point) in postgres data type?

I'm using:

PostgreSQL 8.4
Postgis-1.5.3 and I find that there is 2 types of geography in there:

'$libdir/postgis-1.5', 'geography_from_geometry' (from function geography(geometry))
'$libdir/postgis-1.5', 'geography_enforce_typmod' (from function geography(geography, integer, boolean))

GeoServer-2.1.3


texts:
I've created a network topology using QGIS, saved it as a shp, imported it into PostGIS using shp2pgsql and I'm now trying to output some of the geometry using mapserver.

Unfortunately mapserver is kicking up the following error:

  msDrawMap(): Image handling error. Failed to draw layer named 'ROUTE'.
  msPostGISLayerWhichShapes(): Query error. Error (ERROR: Operation on
  two geometries with different SRIDs ) executing query: select
  encode(AsBinary(force_collection(force_2d("the_geom")),'NDR'),'hex')
  as geom,"gid" from (SELECT buildingroute.the_geom,buildingroute.gid
  FROM routetest INNER JOIN buildingroute ON
  routetest.toid=buildingroute.gid ORDER BY routetest.id ASC) AS foo
  where the_geom && GeomFromText('POLYGON((-22000000 -11000000,-22000000
  11000000,22000000 11000000,22000000 -11000000,-22000000
  -11000000))',32709)

Now I don't know a lot about spatial referencing, but I've copied the settings from QGIS and I still get those errors. Here's those settings in QGIS:

The mapserver file I created can be found at http://pastie.org/2202322

If anyone can give me advice on what I need to change in the mapserver file I'd appreciate it.. I expect I've just got the projection settings wrong..

Thanks!


texts: Any performance difference between Python window and Python script tool?

texts:
Is there a difference in performance between running a Python script in the Python window versus running it as a Python script tool from the toolbox in ArcGIS Desktop 10.0?

Sometimes I need to use ArcPy and have ArcGIS running. In my recent experience I can run the same script via toolbox and in the python window and the toolbox seems to be a bit faster.


texts: Generate Near Table causes ArcMap shutdown

texts:
When I try to find the distances from xy coordinates to linear landforms such as roads and rivers, using Generate Near Table, it keeps shutting arcmap down and sending error reports, without telling me what I have been doing wrong?
Does anyone know a way around this? Or what I am doing wrong? I'm pretty certain that generate near table is the only way to get the data that I need.


texts: Where are descriptions of attributes in datafiles on naturalearthdata.com?

texts:
I have downloaded some datafiles from naturalearthdata.com site.

For example, ne_50m_admin_0_countries feature source contains dozenz of attributes with names like scalerank, level, mapcolor7, mapcolor8, mapcolor13 and so on.

It is possible to guess the meaning of some of them but others are mysterious.

Is there any description for these attribute names?


texts: Drawing 3d polygon element on Arcobject axSceneControl?

texts:
I have standalone application (ArcObjects C# 10.1) uses axSceneControl. How to display on the scene 3d polygon elements?
I have implemented drawing 3D line elements using ISimpleLIne3DSymbol. I tried to use PolygonElement to add to AxSceneControl, but the element has no elevation.

long pColor = Vsur.GIS.Element.ColorConvert.ConvertColor(System.Drawing.Color.Blue);
IElement pElement = ElementOperate.CreatePolygonElement(geo, pColor, esriSimpleFillStyle.esriSFSDiagonalCross, pColor, esriSimpleLineStyle.esriSLSNull, 0.5, false, "OPERA");
IGraphicsContainer3D pGCon3D = pSceneControl.Scene.BasicGraphicsLayer as IGraphicsContainer3D;
pGCon3D.AddElement(pElement);


texts:
I'm working in the SRID 4326 system and need to compute the area of rectilinear polygons.  Our geometry (geom) column is in WKB format

When I execute the ST_Area function it appears my results are returned in squared degrees when I execute a call as follows:

public.ST_Area(geom)/1000000 from production.schemaA where item = '123456'

Then if I cast the geom column to geography the result returned was 33.97 km2

This is still quite a bit off from what I get checking the area of the SAME polygon in Arc Map and QGIS Desktop (Field calculator method) about 43.37 km2

I have investigated using the ST_Transform function to transform out of 4326 but this doesn't work as you have to know which UTM/WGS SRID to use in advance based on the (UTM/ZONE) zone you are working in. We'll be hitting our global database of polygons all over the world had identifying the specific zone is cumbersome.

Example:
public.ST_Area(geom::geography)/1000000 from production.schemaA where item = '123456'

There are threads that are close to this issue but not exactly the same issue of have no solution.


texts: Batch clipping Shapefiles using Argparse

texts: Is there a way to batch clip shapefiles in ArcPy using an argument-Argparse (which calls a shapefile or group of shapefiles) specified by an outside user in command line?


texts: Converting raster layer to vector using ArcGIS Desktop?

texts: I have a raster layer (processed from landsat thermal band). But this data does not contain attribute table. I want to convert this into a vector layer of polygon features; holding an attribute table of temperature data. So I need to convert this to contain an attribute table. I am using ArcGIS software.


texts: Ways to change ArcMap 10 editing tool back to ArcMap 9.3 version

texts:
I'm trying to get used to using ArcGIS 10, but in the meantime I still need to actually get work done.  I stumbled across a way to change the editor tool from its 10 version back to the 9.3 version (with the drop-down task bar etc).  Now I can't find out ow to do it anywhere, -1 internet for me.

Can anyone can give me the steps to change the editor toolbar back to the 9.3 version?


texts: Giving color ramp for our data based on values using cartodb api

texts:
I have average values in the data stored in an array. Based on these average values, does CartoDB have functionality to fill the polygons in different color?Like if I give minimum and maximum values of average in a function it should fill colors automatically.

I am sure it can be done using the CartoDB editor. Can this be done using the CartoDB APIs?


texts: Recovering accidental edits through ArcGIS Pro Version Changes window?

texts:
We have a situation with our ArcSDE where features were accidentally deleted inside a version and we would like to recover those changes. We cannot undo edits because they were saved in a previous session, we have reconciled the version once, and we do have backups of the data if that's what it comes down to.

It seems possible to recover accidental edits through ArcGIS Pro's Version Changes window, but is there a particular workflow that needs to be followed to make replacements stick?

We know this was not possible in ArcMap, but playing around in Pro suggests that the same functionality we've had for resolving conflicts post-reconcile has been carried over to Version Changes:

In this particular situation, we want to replace "Current" (the accidental change that happened in the version) with target (what's shown in DEFAULT). However, when we tell the application to "Replace with Target Version", nothing resolves.

My extensive research has turned up very little so far.

We're running Pro 2.2.4 with Oracle SDE, and our DBA tells me he compresses states twice a week.


texts: How can I deploy a custom WPS on deegree without using the UI?


texts: How do I connect to geoserver on my system using java code?

texts:
I am writing the following code...
But i am getting an error as...

java.io.FileNotFoundException: http:// localhost:8080/geoserver/
    at sun.net.www.protocol.http.HttpURLConnection.getInputStream(HttpURLConnection.java:1613)
    at Try.main(Try.java:25)

Am i doing it right? I am new at this geoserver as well as client server coding...

import java.io.*;
import java.net.*;

public class Try{

    public static void main(String args[]) throws Exception{

        URL hp= new URL ("http://localhost:8080/geoserver");

        try {

            URL myURL = new URL("http://localhost:8080/geoserver/");

            URLConnection myURLConnection = myURL.openConnection();

            myURLConnection.connect();

            BufferedReader in = new BufferedReader(new InputStreamReader (myURLConnection.getInputStream()));

            String inputLine;

            while((inputLine = in.readLine()) != null )
            {
                System.out.println(inputLine);

            }

            in.close();

        }
        catch (Exception e){
            System.out.println("Into CATCH");
            e.printStackTrace();
        }
    }
}


texts: Assigning a new coordinate system and projection to a shapefile in R

texts:
I have a shapefile with no spatial reference nor projection, and a GEOTIFF that has the spatial reference and projection data. The shapefile was created from a segmentation analysis of the GEOTIFF, therefore the extent of the shapefile matches that of the GEOTIFF, and should be overlayed on the GEOTIFF. The two files can be found here.

I have tried the raster, sp and RGDAL packages, among others, to try to project/reproject the shapefile. Due to the size of the shapefile, it takes over 2 days to rasterize from the raster package. I have tried to convert a smaller shapefile (samp.shp) to a raster, then copy the projection of the GEOTIFF (geotiff) to the new sampled raster.

samp.shp.raster<-raster(samp.shp,res=r)
shp.rasPRO = projectRaster(samp.shp.raster, geotiff)

This seems to work on the sample, but I cannot get to reproduce it for my data. I think if I can convert the shapefile to a grid, then I can copy the X,Y coordinates of the geotiff to that grid, but I am not sure if this would work. I have >200 individual GEOTIFFs and corresponding shapefiles and I need to batch process. I am searching for the most efficient method of assigning the spatial reference and projection data to the shapefiles.

Any suggestions would be welcomed.


texts: Changing pixel values in Google Earth Engine?

texts: What is the best method to change the negative values (or values greater than some number) in every band of a multiband image to zero?


texts: Inconsistent pixel values shown on openlayers map

texts:
I'm getting weird behavior with different zoom levels of my WMS layer. It either omits data completely, or changes pixel data and I don't see where the error/misconfiguration could reside.

Openlayers layer:

var chl = new ol.layer.Tile({
          opacity: .5,
          source: new ol.source.TileWMS({
            attributions: 'Chlorophyll Satelite Data',
            url: "/cgi-bin/mapserv.exe?map=\\ms4w\\apps\\tinyows\\myDemo\\WMS.map",
            params: {'LAYERS': 'CHLLLL', tiled : true },
            serverType: 'mapserver',
            wrapX: false,
          })
        });

Map File:

LAYER
  NAME "CHLLLL"
  DATA ".\data\CHL\chlor_a_tiled.tif"
  EXTENT -180 -90 180 90
  TYPE RASTER
  PROJECTION
    "proj=latlong"
    "ellps=WGS84"
    "datum=WGS84"
  END
  STATUS ON

   METADATA
    wms_title "Chl"
    wms_abstract "Chl"
    wms_srs "EPSG:4326"
    wms_include_items "all"
    wms_enable_request "*"
   END

        CLASS
          EXPRESSION ([pixel] < .001)
          STYLE
            BACKGROUNDCOLOR 0 0 0
          END
        END
        CLASS
          EXPRESSION ([pixel] > .01 AND [pixel] < .1)
          COLOR 237 248 233 #GREEN
        END
        CLASS
          EXPRESSION ([pixel] > .1 AND [pixel] < 1)
          COLOR 186 228 179 #GREEN
        END
        CLASS
          EXPRESSION ([pixel] > 1 AND [pixel] < 10)
          COLOR 116 196 118 #GREEN
        END
        CLASS
          EXPRESSION ([pixel] > 10 AND [pixel] < 30)
          COLOR 49 163 84
        END
        CLASS
          EXPRESSION ([pixel] >= 30)
          COLOR 0 109 44 #GREEN
        END
END

And how I prepared my tif from nc:

gdal_translate -of Gtiff -a_srs EPSG:4326 HDF5:"V20170012017008.L3m_8D_SNPP_CHL_chlor_a_4km.nc"://chlor_a chlor_a.tif
#Tried seeing if this made a difference
gdal_translate -of Gtiff -co "TILED=YES" chlor_a.tif chlor_a_tiled.tif

Map declaration:

const map = new ol.Map({
  layers: [
  new ol.layer.Tile({
    source: new ol.source.TileArcGISRest({
      url: 'http://services.arcgisonline.com/arcgis/rest/services/Ocean/World_Ocean_Base/MapServer',
      wrapX: false
    })
  }),
  moorings,
  chl,
  cones,
  finWhaleCoverage,
  //sst,
  //ais,
  ],
  target: 'map',
  view: new ol.View({
    projection: 'EPSG:4326',
    center: [-78,34],
    minZoom: 4,
    maxZoom: 8,
    zoom: 5.5,
    extent: [-100,-60,30,60],

  })
});


texts: Openlayers support for vector tiles, pbf split and crop

texts:
I have an Ionic 2 app with an embedded map, it needs to work 100% offline so I chose to use openlayers 3 with local z/x/y tiles in png format. Now, the problem with that is the size of the map.
I'm thinking about switching to a vector tile, but it needs to be offline as well. I know that OL3 supports pbf vector tiles servers

var layer = new ol.layer.VectorTile({
source: new ol.source.VectorTile({
    attributions: '© OpenMapTiles ' +
    '© ' +
    'OpenStreetMap contributors',
    format: new ol.format.MVT(),
    tileGrid: tilegrid,
    tilePixelRatio: 8,
    url: 'https://free-0.tilehosting.com/data/v3/{z}/{x}/{y}.pbf?key=tXiQqN3lIgskyDErJCeY'
})
});

But I never found an example of tile in a local pbf (or any other vectorial format), nor a way to split and crop a big pbf to only include the interesting area.

Are local pbf supported?

Can I use something like
url: 'assets/tiles/{z}/{x}/{y}.pbf'
to include a local vector file?

How do I split and crop to get the {z}/{x}/{y} folder structure only for the relevant area?


texts: How to trigger registered events in Openlayers?

texts:
I am building a web application where I need to compare 2 different maps of the same area, I need for the events on one map to automatically be triggered on the other map.
I have so far: map2.events.triggerEvent(evt.type,evt.px);

these guys have done a map compare: http://tools.geofabrik.de/mc/?lon=29.53125&lat=22.26877&zoom=2

but I have no idea bow they achieved this event handling..

please help, I need this functionality As soon as possible, its been a week am fighting with this..I will appreciate all replies.


texts: How to automatically generate a route between locations?

texts:
I need to create a route between 70-80 different locations or points in a map, I would like something a la google map, I select my origin and destination and it automatically generates a route between the points

I can't use google maps because it has a limit of 10 destinations per route. And so far I haven't been able to do it with QGIS. Is there any plug-in or way to do what I need?


texts: How can I generate a high resolution rendering of the globe?

texts:
I am building an ANSI E poster size map collection and would like to add a nice rendering of the globe to the maps indicating the approximate position of the region on the globe.

The problem is that I need the globe rendering to be high quality, high resolution image. I think minimum 300 or 400DPI @ about 12" x 12" inches so that I can work with it in CorelDRAW where I will be creating the final product.

I have Google Earth (Free), NASA World Wind and ArcGLobe but neither software is capable of producing a satisfactory result. All look OK but the resolution of the image is not a press quality.

So far Google Earth seems best (although terrain relief would be better then the spotty imagery), followed by NASA World Wind and legging behind is ArcGlobe. (Given the price tag, one would expect ArcGlobe to deliver but I am fully aware of the ArcScene exporter problems so I am not very surprised.)

If someone could recommend a virtual globe (free) program that I can use to do this that would be great. Other options are welcomed as well; the globe does not have to be dynamic but it has to focus on Canada, centered on Ontario.

Here are the outputs from the above mentioned programs:

Google Earth

NASA World Wind

ESRI ArcGlobe 10


texts: SPOT6 image does not display correctly in arcpad

texts: I am using  ArcPad Data Manager from Arcmap 10.3 to export SPOT6 image to .apm file for Arcpad, and using it for field work, I exported it as background TIFF, the projection I used is :WGS_84_UTM_zone_23N.  Everything was fine, but the image display somehow changed, I thought it because of SPOT6 has 4 bands, I composite them with first three bands, I also exported the image as GRID format, but none of them works. You can see clear difference between two images(below).Can anyone tell me how to figure this out?


texts: Why is ArcGIS Server so slow in creating mapbook PNG

texts:
We have 3 MXDs for creating 3 different mapbooks. Our customer want a daily export to PNG. For this, I use a Python-Script, which works fine for me on my Desktop. But I need it on ArcGIS Server. There it works also, but is really slow. In fact, it takes up to 2.5 times on Server than on Desktop.
Two of the MXDs have to export 173 PNG, the third exports only 32.

Times on Desktop: One of the 'big' MXD (with 173 PNG) needs ca. 5 hours ( ~106 seconds per PNG), but on the Server it needs ca. 12 hours ( ~255 seconds per PNG).

I tried a SDE as datasource (which is the way i normally want to go) and a fGDB. Also Labelling on and off. But there where no specific difference in processing speed.

Server is an ArcGIS Server 10.2.1.

Can anyone give tips for a better processing time or have the same issue? And please, know comments like these mapbooks are too big for a daily export ;)

Another question, if someone have experience with that: How much could multithreading give here?


texts: Asp.net or MsSql Create bbox of a point proportional to image height and width

texts:
Short question is how do I create bbox dimensions that is proportional to image height and width?

I'm working on a asp.net webservice to generate a map using geoserver. The parameters I have are Point, Width and Height of Image, and zoom level.

My need to calculate the bbox size and position based on the point and the height and width of the parameters.

    public void MapCreate( string Point, int Width, int Height, int ZoomLevel )
    {

        // Calculate bbox using sql or c#
        string bbox = "-12942776.251415756,5411282.980420768,-12942623.377359185,5411435.8544773385";

        // Url for GeoServer
        string url = string.Format( "/wms?service=WMS&version=1.1.0&request=GetMap&layers=vvv:ccc&styles=&bbox={0}&width=432&height=468&srs=EPSG:3857&format=image/png", bbox, Width, Height );
    }

I was thinking it might be easier to do this in mssql but I have no idea on how to pursue on this.

   -- SP Parameters
   DECLARE @ImgWidth INT;
   DECLARE @ImgHeight INT;
   DECLARE @Zoom INT;
   SET @ImgWidth = 400;
   SET @ImgHeight = 600;
   SET @Zoom = 10;

   DECLARE @p geography;
   SET @p = geography::Point(40, -40, 4326);

   SELECT @p
   UNION ALL
   SELECT geography::STPolyFromText('POLYGON((-122.358 47.653, -122.348 47.649, -122.348 47.658, -122.358 47.658, -122.358 47.653))', 4326); -- Calculated bbox


texts:
I've been searching for a way to generate server side static maps using Python (or any other languages, doesn't matter).

I can create the maps I want with QGIS (and they look really great) but I don't know how to automate the process.

What I want is;

Creating a map style in QGIS and export it to some kind of XML (there is a plugin named quantumnik but it is not available any more I guess)
Running a query over PostGIS in my Python script
Creating a shapefile or only a PNG using the dataset returned from query and style XML defined above.
Repeating the whole process every night with different datasets.

I've played around with Mapnik but couldn't succeed in generating a heatmap yet (https://github.com/stellaeof/mapnik-rasterizers)

I've read about QGIS server, GeoServer, MapServer. But now I am confused since there are a lot of tools around.

Is Mapnik the right direction or are there any better tools for that kind of process?
Or is there a way to use QGIS programmatically on a server?

EDIT:

After trying a lot of different things, I've decided to use GeoServer. I've installed it using using docker container image: https://github.com/kartoza/docker-geoserver
It can generate png, gif, jpeg, and also can serve WMS layers that can be integrated with Leafletjs pretty easily.

Only problem is, I am trying to design my maps with QGIS first and then export them to GeoServer using "Geoserver Explorer" plugin and it doesn't work well for labels and some other features. There is no easy way to generate SLD in QGIS I guess. I keep fixing SLD configurations by hand after exporting them to GeoServer. But this is still the best solution I guess.

It would still be great if you can suggest better plugins or ways to integrate QGIS and GeoServer or any other SLD editors.


texts:
I am using proj, in particular pyproj, to convert from UTM coordinates to the (x,y) coordinates used, for example, by Google Maps. There is a parameter called zone which I could not understand how I can determine. I have seen it varying across examples I looked up. What is the rule for specifying the zone?

By the way, the reason I need this transformation is because the points my shp are in this format once I use pyshp to access them.


texts:
I know Solaris is on the way out with 10.1

I also know that v10 is only supported on 32bit flavours of RedHat & SUSE.

What I do want to know, has anyone had any luck installing it on either Ubuntu or CentOS?

Google only told me of this blog which had revolved around Engine 9.2.

--EDIT--
Looks like Ubuntu support for may arrive at 10.1.  Ease of use and lower costs on AWS look to be the reasons (over CentOS)


texts: Zooming to selected features from multiple layers using QGIS

texts: Is it possible to zoom to features selected from multiple layers in QGIS?


texts: Distributing Python script tools using methods not in Sharing Tools section of ArcGIS Help?

texts:
I want to be able to distribute a custom script and it seems that I cannot create a .tbx file, so I looked into the other methods, but they are all only for Python modules it seems.

What I want to be able to do is have the user still be able to pull up the dialog box prompting for input when they click on the tool, how do I go about packaging the external libraries that I used for this script and distribute it such that the user can simply import the tool into his toolbox and click on it and it will pull up the dialog box for user input?


texts: How to convert .csv with WKB to shapefile (.shp) using open source solutions?

texts:
I have a .csv file with hundreds of rows pulled out via PostGIS.

The geospatial information is in "the_geom" field which seems to be in WKB-format.

How can I convert this to shapefile using simple free tools?


texts: Circular focal means with raster and NAs

texts:
I am trying to understand circular focal means using rasters having NA values. I made a couple of scenarios to better understand how to combine weight matrices and functions to summarize values. I am using a circular weight matrix considering the focal cell and four neighboring cells.

library(raster)

r<-raster(extent(0,1,0,1),ncol=5,nrow=5,vals=1:25)
wm<-focalWeight(r,d=res(r)[1]*1.001,type="circle") # original weight matrix
w<-wm
w[w>0]<-1

w
     [,1] [,2] [,3]
[1,]    0    1    0
[2,]    1    1    1
[3,]    0    1    0

Using a first example without NA values and if we focus on the middle value of the raster, (A) I'm not able to correctly compute it using mean and weights of 0/1. It seems that 0 values are included in the calculation as the denominator in the mean calculation is the number of values in the weight matrix (9). The documentation warns about the fact that using weights different from 1 in combination with na.rm = TRUE may lead to problems in calculations, but in this case I would have expected to get the correct value. (B) I can compute the correct value by using instead weights of 0/0.2 and the sum function.

par(mar=c(1,1,5,1),mfrow=c(2,2))
plot(r,axes=FALSE,legend=FALSE,main="Initial raster")
text(coordinates(r),labels=getValues(r))

### A
w<-wm
w[w>0]<-1
r2<-focal(r,w=w,fun=mean)
plot(r2,axes=FALSE,legend=FALSE,main="A\nMean\nweights = 0/1")
text(coordinates(r2),labels=round(getValues(r2),2))
mtext(side=3,line=-1,"(8+12+13+14+18) / 9 = 7.22",cex=0.7)

### B
w<-wm
r2<-focal(r,w=w,fun=sum)
plot(r2,axes=FALSE,legend=FALSE,main="B\nSum\nweights = 0/0.2")
text(coordinates(r2),labels=round(getValues(r2),2))
mtext(side=3,line=-1,"(8+12+13+14+18) * 0.2 = 13",cex=0.7)

In this second example, I turn the middle value to NA and try to compute the missing value using the same focal area. As in the first example, the mean method (C) does not work and the denominator is now 8 since the missing value was removed. The sum method (D) does not work anymore, because the weights are not adjusted for the removed NA value. A solution (E) could be to use a 0/1 matrix and my own mean function to remove all 0 weight values before the mean calculation.

r[13]<-NA

par(mar=c(1,1,5,1),mfrow=c(2,2))
plot(r,axes=FALSE,legend=FALSE,main="Initial raster with NA")
text(coordinates(r),labels=getValues(r))

### C
w<-wm
w[w>0]<-1
r2<-focal(r,w=w,fun=mean,na.rm=TRUE)
plot(r2,axes=FALSE,legend=FALSE,main="C\nMean with na.rm = TRUE\nweights = 0/1")
text(coordinates(r2),labels=round(getValues(r2),2))
mtext(side=3,line=-1,"(8+12+14+18) / 8 = 6.5",cex=0.7)

### D
w<-wm
r2<-focal(r,w=w,fun=sum,na.rm=TRUE)
plot(r2,axes=FALSE,legend=FALSE,main="D\nSum with na.rm = TRUE\nweights = 0/0.2")
text(coordinates(r2),labels=round(getValues(r2),2))
mtext(side=3,line=-1,"(8+12+14+18) * 0.2 = 10.4",cex=0.7)

### E
w<-wm
w[w>0]<-1
r2<-focal(r,w=w,fun=function(i){mean(i[i>0],na.rm=TRUE)})
plot(r2,axes=FALSE,legend=FALSE,main="E\nCustom mean function\nweights = 0/1")
text(coordinates(r2),labels=round(getValues(r2),2))
mtext(side=3,line=-1,"(8+12+14+18) / 4 = 13",cex=0.7)

The problem with this solution is that (1) it relies on my raster not having values of 0 and (2) using my own function greatly increases computation time. The actual matrix I need to run this on is about 13000 x 13000 pixels with lots of NAs.

n<-2000
r<-raster(extent(0,1,0,1),ncol=n,nrow=n,vals=1:(n^2))
w<-wm
w[w>0]<-1
system.time(r2<-focal(r,w=w,fun=mean,na.rm=TRUE))
   user  system elapsed
   0.29    0.00    0.30
system.time(r2<-focal(r,w=w,fun=function(i){mean(i[i>0],na.rm=TRUE)}))
   user  system elapsed
  74.42    0.40   76.28

My question is this: I there a more efficient way than using a custom-made function to correctly compute circular focal means when NA values are present?


texts: Saving results in automap R package for time series data

texts:
I'm planning to interpolate monthly rainfall data from  1980 and 2010. I'm using following code:

require(gstat)
require(sp)
require(automap)
require(rgdal)
require(raster)

for(i in 1:ncol(data@data)) {
  kriging_new <- autoKrige(data[,i],data,grid)
  plot(kriging_new)
 }

I want to save predicted values (kriging estimates), kriging prediction variance, kriging standard deviation for each i and save it as a img file. Can anyone help me to set up this within the loop.


texts: how to list table rows of a layer in a listbox?

texts:
I'm new in using c# and Arcengine .
 i would like to list table rows of a feature in a Listbox like select by attribute in ArcGIS
such that when i click a button, table rows list in a listbox . i use two combobox to select layers and fields . how can i write this code in C# using Arc engine? i don't know how do i write loop for this purpose. i only know i want to run the application in ArcGIS then i use Ihookhelper and Itable.


texts: GDAL, python: How to create a GDAL dataset object from a numpy array

texts:
I am working on a program that uses both OpenCV and GDAL to perform different tasks. OpenCV stores its images as numpy arrays. I would like to be able to go directly from this numpy array to a dataset object without having to save the image to the file system and then reopen it again in GDAL. Right now, I have to do this for each image:

import cv2
import numpy as np
from osgeo import gdal

img = cv2.imread("image.png", flags=-1)
# modify the image
cv2.imwrite("output.png", img)
ds = gdal.Open("output.png")
# do stuff with GDAL

The issue is that the cv2.imwrite step is very time consuming, especially since I'm doing it for a lot of images.

My question is this: is there a way to get a dataset from a numpy array directly?

Something like this:

import cv2
import numpy as np
from osgeo import gdal

img = cv2.imread("image.png", flags=-1)
# modify the image
ds = gdal.datasetFromArray(img)
# do stuff with GDAL


texts: Access attribute table of in-memory raster

texts:
I'd like to read from the attribute table of an intermediate raster generated by arcpy.sa.Combine:

combined = arcpy.sa.Combine(['raster1', 'raster2'])
arcpy.BuildRasterAttributeTable_management(combined)
data = [row for row in arcpy.da.SearchCursor(combined, ["VALUE", "COUNT"])]

The above code raises:
     RuntimeError: 'in_table' is not a table or a featureclass

Is there a way to access the RAT without saving the raster?


texts: Control rendering/layer order of outputs from QGIS 3 Processing Modeler

texts:
I'm using QGIS 3.0 model processor. I am trying to control the order in which the outputs render. Any ideas? I have tried to use parent algorithms but it seems to not make a difference. Is this a bug or an incorrect use of the parent algorithm function. Pictured are the outputs from my model. I want to have all the "A" files automatically draw on top of the others when the output from the model. ]

I've started trying to work with the code given by Kadir below and integrating it into the below hyperlink he provided. I have able to load in and work the example code. Now I am working on integrating Kadir's code into the model. I've included a snippet of the model below as well.
https://anitagraser.com/2018/03/25/processing-script-template-for-qgis3/


texts: "Geometry" Problem while creating Postgis database on Ubuntu

texts:
I am trying to build pgrouting on ubuntu. Therefore, I now firstly try to create a Postgis database on ubuntu.

I am following the instructions on this site:

Installing PostGIS 1.5 on PostgreSQL 8.4 on Ubuntu

I have installed postgresql and postgis.

I now try to run the necessary sql queries:
psql -d my_db -f /usr/share/postgresql/8.4/contrib/postgis-1.5/postgis.sql

However, I got an error as shown in the screenshot below:

http://postimage.org/image/ggk7xg17p/

It looks like there is a problem with the query as
-it firstly aborts the transactions continuously and
-then says type "geometry", "box2d", "geometry" cannot be found

I was expecting these variables to occur as a result of this query. Why am I getting such an error while running the query?

Any idea would be more than welcome.

Thanks!

I checked the first error I get on this action:

ERROR Could not access the file "$libdir/postgis-1.5": No such file or directory

I have been checking on Google but not so useful yet. Any idea to fix it?


texts: How to set GDAL/cs2cs/proj to use NTv2 grid shift file when using EPSG codes?

texts:
I am using the GDAL API for Python for a script that takes only EPSG codes for inputs to specify source and target coordinate systems.

While the cs2cs command seems to handle quite well transformation between NAD27 (EPSG:4267) and NAD83 (EPSG:4326) coordinates in North America, I would like to have similar capabilities for AGD66 (EPSG:4202) to GDA94 (EPSG:4283) and AGD84 (EPSG:4203) to GDA94 (EPSG:4283) in Australia.

When using the -v option with cs2cs to look at the proj string being used i get the following:

    >cs2cs -v +init=epsg:4267
    ># ---- From Coordinate System ----
    >#Lat/long (Geodetic alias)
    >#
    ># +init=epsg:4267 +proj=longlat +datum=NAD27 +no_defs +ellps=clrk66
    ># +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat
    ># ---- To Coordinate System ----
    >#Lat/long (Geodetic alias)
    >#
    ># +proj=latlong +datum=NAD27 +ellps=clrk66
    ># +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat

The nadgrids parameters tells cs2cs to use the best fitting datum grid shift file from conus, alaska, ntv2_0.gsb, and ntv1_can.dat.

For AGD66 (EPSG:4202) cs2cs returns:

    >cs2cs -v +init=epsg:4202
    ># ---- From Coordinate System ----
    >#Lat/long (Geodetic alias)
    >#
    ># +init=epsg:4202 +proj=longlat +ellps=aust_SA
    ># +towgs84=-117.808,-51.536,137.784,0.303,0.446,0.234,-0.29 +no_defs
    ># ---- To Coordinate System ----
    >#Lat/long (Geodetic alias)
    >#
    ># +proj=latlong +ellps=aust_SA
    ># +towgs84=-117.808,-51.536,137.784,0.303,0.446,0.234,-0.29

There is no nadgrids parameter but a towgs84 instead.
A grid shift file from the ICSM; "A66 National (13.09.01).gsb" contains the shift between AGD66 and GDA94.

How can I make GDAL/cs2cs/proj identify this file when using EPSG codes?

I tried to reinstall GDAL/cs2cs/proj using the osgeo4w installer with A66 National (13.09.01).gsb within the directory containing the other datum shift files prior to the reinstallation but without success.


texts: Deleting overlapping polygon features from one of two layers using QGIS?

texts:
I have two overlapping layers of polygons, and I need to somehow erase the areas the overlaps. -This way i should end up with two layers of non overlapping polygons. Sounds pretty simple but I just can't get it right with any of the tools in Qgis!?
Should i use "clip" or "cut shape by polygons", or perhaps the "difference" tool and the "clip"?

Whatever I do the product isn't right.


texts: Sorting (ORDER BY) rows from Shapefile using Python?

texts:
I'm writing a Python script where I have to read a Shapefile, sort its rows/records by a specific field, and then write a new shapefile with that new order.

Does anyone how can I do it?

I'm using Fiona to open the file but I can't find out how I can do an so easy and common operation as an ORDER BY.

EDIT

The solution proposed by @Jose works perfectly right. But geopandas doesn't create the .prj file that contains the projection description and this is a problem if you need it. A possible solution could be just copy the original file with the new name:

import os
import geopandas as gpd

def main(file, field):
    output_file_name = os.path.splitext(file)[0] + '_sorted_by_' + field + os.path.splitext(file)[1]

    shape = gpd.read_file(file)
    shape_sorted = shape.iloc[shape[field].sort_values().index.values]
    shape_sorted.to_file(driver='ESRI Shapefile', filename=output_file_name)

    # Copy prj file
    prj_file = os.path.splitext(file)[0] + '.prj'

    if os.path.isfile(prj_file):
        from shutil import copyfile
        output_prj_file_name = os.path.splitext(file)[0] + '_sorted_by_' + field + '.prj'
        copyfile(prj_file, output_prj_file_name)

Finally I decided to use a ogr2ogr command with -sql argument and execute it in Python using os.system(command) utility.

import os

def main(file, field):
    output_file_name = os.path.splitext(file)[0] + '_sorted_by_' + field + os.path.splitext(file)[1]
    file_object = os.path.splitext(os.path.basename(file))[0]

    command = 'ogr2ogr -sql "SELECT * FROM ' + file_object + ' ORDER BY ' + field + '" ' + output_file_name + ' ' + file

    os.system(command)

if __name__ == '__main__':
   main()


texts: How to compute and fill One-way field in a street shapefile according spacial position

texts:
I think it is something complicate, but I hope someone could bring up any idea about.

I'm using ArcGIS10, and I have a "One-way roads" shapefile with a One-way field with Y/N. But I need to create a new Field for this layer containing F or T if:

F : The final point of the line is Northern and Western than the origin of the line.
T : The final point of the line is Southern and Eastern than the origin of the line.

An example,

Another Example:

Thanks for your time,


texts: QGIS freezes when trying to view layer properties

texts: As of recently whenever I try to edit the properties of a vector layer the entire program just freezes and I'm forced to end the process using task manager. I use Windows 10. Any reasons as to why this might be happening?


texts:
I want to visualize direction of river's course by an arrow placed next to the river's name label.
I have an idea to place the blue arrow as a text in front of my label or at the end of my label, but the problem is with the qgis label engine. Depending on vector rotation label puts along or across the vector direction. Can I check somehow if label is placed along or across the vector direction? I would prefer the solution by Expression or maybe you have any better ideas how to visualize that...
Expected result below: (I marked the river's course (which equals vector direction) with the red arrow):


texts: Binding multiple Leaflet image overlays as one

texts:
I am planning to put at least ten single image overlays, spread across the country. I have got two so far:

var priorycourt = 'Priory_court.PNG',
imageBounds = [[52.399892, -1.531154], [52.400837, -1.52945]];
L.imageOverlay(priorycourt, imageBounds).setOpacity(0.75).addTo(map);

var ipswich3a = 'A3Ipswich.svg',
imageBounds = [[52.0303, 1.19047], [52.0326, 1.19624]];
L.imageOverlay(ipswich3a, imageBounds).setOpacity(0.75).addTo(map);

Which won't look good, when it is 70 instead of 2.

I tried to put them in a different js file img.js,
although only very first is working, when in my main script remain

 L.imageOverlay(priorycourt, imageBounds).setOpacity(0.75).addTo(map);

Is it possible to create something analogue to L.LayerGroup[priorycourt,ipswich3a]... in order to combine all these image overlays together in one piece?


texts:
I am using GDAL 1.11.3 to merge about 9000 tif of satellite images (for a total of 60 GB) to create a huge image, representing the whole Mediterranean, which I will then need to transform into XYZ tiles for my TileServer. These tiffs have different projections, so I transform them all into my projection:
gdalwarp -t_srs EPSG:4326 -r average -dstalpha %mypath%!infile! %mypath%!outfile!

When I launch the merge, the application freezes after a few hours. Maybe the 60GB of final image needs to be loaded into RAM, running out of it?
python3 -m gdal_merge -ot Byte -of GTiff -o c:\merged.tif --optfile c:\tiffsList.txt

Is this the right way to process these tiffs to build a mosaic?


texts: Debugging arcpy.da.UpdateCursor where clause which gives RuntimeError: Cannot find field?

texts:
I'm trying to use UpdateCursor in arcpy, but the where clause I'm using seems to be causing trouble. However, when I open up the table in ArcMap and use Select By Attributes with the exact same clause that causes an error in the Python script, it works fine.

The name of a feature class and a dictionary containing strings is passed into a function with the following code:

rows = arcpy.da.UpdateCursor(featureClass, "{0} = {1}".format("PropCode", "'"+hotelDict["hotelId"]+"'"))
for row in rows:
    print "stuff" #I just wanted to see if it was making it here

I get the following error: RuntimeError: Cannot find field 'PropCode = 'FR066''

More info that may or may not be helpful: I'm using ArcGIS 10.1, the field that is being searched in contains text, and using arcpy.AddFiedlDelimiters(featureClass, "PropCode") instead of just "PropCode" did not fix the problem.

There are no results from Googling "RuntimeError: Cannot find field", and PropCode = 'FR066' is correct when using Select By Attributes in ArcMap, so I'm kind of stumped.


texts: How to Geocode 300,000 addresses on the fly?

texts:
I have a database which has 300,000 address, which are to be shown on map. I know if I geo-code all the address it will be too expensive for me. So I was wondering if it is possible to geo-code the address on the fly/ realtime, when a user would select a address (a property address) it would search through the database and then geo-code the address and then map it with other attributes.

It would be really nice if you could share a code, concept, or anything. By the way my backend is on mysql supported by Joomla.


texts: split lines into segments without creating new layer in QGIS

texts: How to split line objects into segments/multiple objects based on either polygons or lines without creating a new layer? The problem is clarified in the picture below. I need the green line to be segmented according to the pink lines or green polygons.


texts: How to fix "TabError: inconsistent use of tabs and spaces in indentation" in QGIS 3.x experimental plugin?

texts:
I have reviewed QGIS 2.18 plugin to port in QGIS 3.x version, but when I run that QGIS 3.x plugin I get this error message:

  File "C:/Users/kapil/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\polygon_divider\Polygon_Divider.py", line 1331
    layers = [tree_layer.layer()
                                ^
TabError: inconsistent use of tabs and spaces in indentation

So, how to solve or fix that error?

My code link is below:
https://github.com/kapildevadk/Polygon-Divider/blob/master/Polygon_Divider.py
Error lines are:

def run(self):

"""

* Run method that performs all the real work

"""

# Fetch the currently loaded layers

layers = [tree_layer.layer()

for tree_layer in QgsProject.instance().layerTreeRoot().findLayers()]

And  extra lines are;

def unload(self):

    """Removes the plugin menu item and icon from QGIS GUI."""

    for action in self.actions:

        self.iface.removePluginMenu(

            self.tr(u'&Polygon Divider'),

            action)

        self.iface.removeToolBarIcon(action)

    # remove the toolbar

    del self.toolbar

def select_output_file(self):

    """

    * JJH: Open file browser

    """

    # get filename from dialog

    filename = QFileDialog.getSaveFileName(self.dlg, "Select output file ","", '*.shp')

    # verify that a name was selected

    if filename != "":

        # clear previous value

        self.dlg.lineEdit_2.clear()

        # make sure that an extension was included

        if filename[-4:] != '.shp':

            filename += '.shp'

        # put the result in the text box on the dialog

        self.dlg.lineEdit_2.setText(filename)

def startWorker(self, inFile, outFilePath, targetArea, absorbFlag, direction):

    """

    * JJH: Run the polygon division in a thread, feed back to progress bar

    """

    worker = ExampleWorker(inFile, outFilePath, targetArea, absorbFlag, direction)

    start_worker(worker, self.iface, 'running the worker')

def run(self):

    """

    * Run method that performs all the real work

    """

    # Fetch the currently loaded layers

    layers = [tree_layer.layer()

    for tree_layer in QgsProject.instance().layerTreeRoot().findLayers()]

    # Clear the contents of the comboBox and lineEdit from previous runs

    self.dlg.comboBox.clear()

    self.dlg.lineEdit.clear()

    # Populate the comboBox with names of all the loaded layers

    self.dlg.comboBox.addItems([layer.name() for layer in layers])

    # populate comboBox_2 with the possible directions

    self.dlg.comboBox_2.clear() # need to clear here or it will add them all again every time the dialog is opened

    self.dlg.comboBox_2.addItems(['left to right', 'right to left', 'bottom to top', 'top to bottom'])

    # Populate the comboBox with names of all the loaded layers

    self.dlg.comboBox_2.addItems([layer.name() for layer in layers])

    # show the dialog

    self.dlg.show()

    # Run the dialog event loop

    result = self.dlg.exec_()


texts: How to use rasterio transform on all bands

texts:
I've imported a 10 band raster tif using rasterio as ds.  I'm attempting to pad the image but keep getting errors with the transform.

print(ds.transform)
pilpadimg = rasterio.pad(ds.read(), transform=ds.transform, pad_width=((0,0),(1,2),(3,4)), mode='reflect')

| 15.00, 0.00, 344117.00|
| 0.00,-15.00, 5082412.00|
| 0.00, 0.00, 1.00|
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-77-9526c3dd4939> in <module>
      1 print(ds.transform)
----> 2 pilpadimg = rasterio.pad(ds.read(), transform=ds.transform.identity(), pad_width=((0,0),(1,2),(3,4)), mode='reflect')

~\AppData\Local\Continuum\anaconda3\envs\pytorch\lib\site-packages\rasterio\__init__.py in pad(array, transform, pad_width, mode, **kwargs)
    279     padded_array = np.pad(array, pad_width, mode, **kwargs)
    280     padded_trans = list(transform)
--> 281     padded_trans[2] -= pad_width * padded_trans[0]
    282     padded_trans[5] -= pad_width * padded_trans[4]
    283     return padded_array, Affine(*padded_trans[:6])

TypeError: can't multiply sequence by non-int of type 'float'

It works great with numpy but there is no transform..

>>> x = np.ones((8,3,3))
>>> np.pad(x, pad_width=((0,0),(1,2),(3,4)), mode='constant', constant_values=0)

pad_width has to be an int and the transform is a float.  I converted the transform list to int and still get the same error

EDIT:
Fixed by simply using numpy to pad instead of rasterio

pilpadimg = np.pad(ds.read(),pad_width=[(0,0),(1,2),(3,4)], mode='reflect')


texts: Methods for importing GeoJSON feature collection into ArcMap?

texts:
I have a GeoJSON file containing a feature collection, with the structure/format validated by http://geojsonlint.com/.

What do people recommend as the best way to import this into ArcMap?  I have about 1000+ files that I need to import and was hoping to find/build a python script to batch import them.

I did try this toolbox:  https://github.com/jasonbot/geojson-madness, but I do get an error informing me that I have "Inconsistent geometry types", despite the the validation website informing me my structure and content is correct.


texts:
I have a geocentric coordinate for the position of an object; I am trying to write a function where by an objects position changes with respect to time - moves across the surface at a constant velocity. p(x,y,t) is the position of the object at time t, and at each time step dt is updated by a constant velocity function v(x, y).

I will try my best to explain how I understand the problem, please correct me if I am wrong.

The transformations cannot be applied directly to the geocentric coordinates, since v(x, y) is describing motion on a flat plane. This origin of this flat plane is located at the geocentric coordinates G(X,Y,Z) and is oriented such that the vector G(0, 0, 0)->G(X, Y, Z) is orthogonal to p(x,,) and p(,y,).

I think the problem can be boiled down to, given a flat space translation (dx,dy) what is the new geocentric coordinate?


texts: Getting incorrect northing/easting in proj4 from converting Lat/long to NAD83 Maine West 1802

texts:
I'm trying to convert latitude/longitude to northing and easting via proj4js.
This is the string im using to project from.
+proj=tmerc +lat_0=42.83333333333334 +lon_0=-70.16666666666667 +k=0.9999666666666667 +x_0=900000.0000000001 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs

And here are the latitude/longitude coordinates:
44.128783,-70.201664

This is what I get back as  a result:
9650143.893240128, -51464079.10841816

But this is not correct. When I try to convert this back to lat/long with this tool:
https://www.earthpoint.us/stateplane.aspx
and the proper zone number 1802 - Maine West (1983)
I get incorrect response. I'm sure I'm using the right zone / project string so why is it shooting out the wrong northing and easting?


texts: Raster overlay without loosing original information

texts: I have been an active ArcGIS user but now moving to QGIS but i miss a tool that allows me to overlay/intersect raster layers e.g soil classes, slope classes, landuse classes but still main this original information in the new single raster layer. So in the new layer i should be able to click any pixel and know the soil, slope and landuse classes. In ArcGIS, combine tool does this by keeping these information in the attributes table through different fields. A QGIS raster doesnt have a attributes table, so whats the alternative solution?


texts: Using ListFields and CalculateField to find max value in variable number of fields to populate another field?

texts:
I need some assistance with the syntax for a CalculateField_management calculation as I have still yet to master all the python syntax rules.

What I am trying to do is find the max value in a variable number of fields in order to populate another field. I am using the ListFields function to discover the desired fields to choose from, but getting that list into the formula is giving me some difficulty.

import arcpy, os, string
LAYERS = arcpy.GetParameterAsText(0)
SLOSHFILEDS = [f.name for f in arcpy.ListFields(LAYERS,"","DOUBLE")
arcpy.CalculateField_management (LAYERS, "MAXSURGE", max(SLOSHFILEDS))

I have tried any number of different string combinations for the max() calc to no avail (not that this particular variation shows that).

Adding/Changing the following to the script doesn't give me the syntax error that I would recieve with the above, but it does give me a "The calculate value is invalid for the row with ObjectID = 0..." x 18,526 (or however many rows are in my table) and then does nothing to the table except populate my MAXSURGE field with 0's.

SLOSHFILEDS = arcpy.ListFields(LAYERS,"","DOUBLE")
fieldNameList = []
for field in SLOSHFILEDS:
        if not field.required:
            fieldNameList.append(field.name)
arcpy.CalculateField_management (LAYERS, "MAXSURGE", max(fieldNameList))

Hard coding the field names into the formula works great, but of course, I will not always have the same number of fields or same field names to work with.


texts: Calculating elipsoidal length of line in PyQGIS

texts:
I am using PyQGIS with QGIS 3.

This code below is what I developed already, what I am doing wrong?

layer = iface.activeLayer()
feature=layer.getFeatures()
for feat in feature:
    feature.setSourceCrs(##What goes here?)
    print (QgsDistanceArea.measureLength(feat))


texts: Constructing Voronoi Polgon and Applying Nearest neigbour in it

texts:
I am learning QGIS 3.XX and currently working on Voronoi Polygons.

I have three layers

Point layer
Building Layer
Population CSV file.

These points layer are basically water tanks and building layer is full of small Polygon in YY area. I need to find which building will take water from the nearest water tanks.

I need to make Voronoi polygon on Point Layer. I tried manually by going through Layer-> Geometry Tools->Voronoi Polygon. And I got the below output.

After doing this, the newly created Voronoi Layer is placed at top of Layer Panel. I need to add this layer at lowest so that all the layers above can be put simultaneously on canvas. By now I am doing manually by placing point layer and building layer above voronoi layer, but how to achieve this using python? My expected output should be

until now my code looks like

from qgis.core import QgsProject
uri='file:///C:/Test/points.csv?
delimiter=,&crs=epsg:4326&xField=Lat&yField=Long'
layer=QgsVectorLayer(uri,'Points','delimitedtext')
QgsProject.instance().addMapLayer(layer)
vectorLyr=QgsVectorLayer("C:/Test/demand/Demand.shp","BuildingLayer","ogr")
vectorLyr.isValid()
QgsProject.instance().addMapLayer(vectorLyr)
uri='file:///C:/Test/data.csv?delimiter=,'
infoLyr=QgsVectorLayer(uri,'Population','delimitedtext')
infoLyr.isValid()
QgsProject.instance().addMapLayer(infoLyr)
csvField='FID_1'
shpField='FID'
joinObject=QgsVectorLayerJoinInfo()
joinObject.setJoinFieldName(csvField)
joinObject.setTargetFieldName(shpField)
joinObject.setJoinLayerId(infoLyr.id())
joinObject.setUsingMemoryCache(True)
joinObject.setJoinLayer(infoLyr)
vectorLyr.addJoin(joinObject)


texts:
In our OpenLayer.js we first call the maximum extent of the map:

"maxExtent: new OpenLayers.Bounds (X1, Y1, X2, Y2),".

After the layer definition follows

"map.addLayers ([lyr1, layer3, layer4]);" and "map.zoomToMaxExtent ();".

In between, the variable for "OpenLayers.Control.GetFeature" is defined. Our problem now is that after our feature query, the map always returns to "MaxExtent".

Does anyone have an idea, how after the feature query the zoom level does not change?

In our OpenLayer.js we first call the maximum extent of the map:

"maxExtent: new OpenLayers.Bounds (X1, Y1, X2, Y2),"

After the layer definition follows "map.addLayers ([lyr1, layer3, layer4]);" and "map.zoomToMaxExtent ();". In between, the variable for "OpenLayers.Control.GetFeature" is defined. Our problem now is that after our feature query, the map always returns to "MaxExtent".

Does anyone have an idea, how after the feature query the zoom level does not change?

After selected feature, the zoom should not change. We only want the "MaxExtent" when loading the map for the first time.


texts: arcpy.SelectLayerByAttribute_management reference a field and a variable

texts:
I have tried every which way but the right way. I can't seem to get the quotes in the right spot. I have a float field called speeds and I want to compare it to a float variable I am calling Min_Speeds.

#arcpy.SelectLayerByAttribute_management(Phase2+"_layer", "NEW_SELECTION", "\"Speeds\" = '{0}'".format(Min_Speeds))
#arcpy.SelectLayerByAttribute_management(Phase2+"_layer", "NEW_SELECTION", float('\Speeds\') >= " float(\Min_Speeds)\)"")
#arcpy.SelectLayerByAttribute_management(Phase2+"_layer", "NEW_SELECTION","float(\Speeds\ >= 4)"
arcpy.SelectLayerByAttribute_management(Phase2+"_layer", "NEW_SELECTION", "Speeds > "[Min_Speeds]")


texts: GRASS: db.in.ogr fails where v.in.ogr succeeds

texts:
This is a follow up on a previous question on the storage of vectors layers imported from PostGIS.

To importing a spatial table from PostGIS I use, sucessfully, the following formulation with v.in.ogr:

> v.in.ogr input="PG:host=pg-server dbname=dsm user=duque004 password=secret" layer="observations.profiles" output=profiles -o
Over-riding projection check
Check if OGR layer <observations.profiles> contains polygons...
ERROR 1: Could not fetch SRS:
 100%
Creating attribute table for layer <observations.profiles>...
Importing 141584 features (OGR layer <observations.profiles>)...
 100%
-----------------------------------------------------
Building topology for vector map <profiles@test>...
Registering primitives...

But if I feed the same parameters to db.in.ogr, it fails to find the table:

> db.in.ogr input="PG:host=pg-server dbname=dsm user=duque004 password=secret" db_table="observations.profiles" output=profiles
ERROR: Unable to open data source <observations.profiles>
ERROR: Input table <observations.profiles> not found or not
       readable

What is going on here?

Edit: these outputs report to GRASS 7.6.0

Update 01-05-2019: This issue is now reported in ticket 3836.


texts: Splitting data from GeoServer layer into multiple layers with OpenLayers

texts:
I have a GeoServer layer that contains all buildings (houses, schools, library, museums) with their locations, data is stored in PostGIS database, and now creating a geoportal for showing that online with some symbols.

Is there a way to split this layer into multiple layers so I can use a built-in function, in OpenLayers, for showing and hiding layers?

I need that because for example if someone wants to see only schools he can turn all other layers off. If this can't be done can someone give me an idea what I should do?

My other idea is to try creating external checkboxes and then use them to show/hide features. I know that I can split them with an SQL view, but there is a lot of data for that. I am using OpenLayers 3.


texts: osmconvert return error "Segmentation fault (core dumped)" when trying to extract a smaller PBF from a bigger one (based on boundaries polygon)

texts:
I have a rather large PBF file and I want to extract a small part out of it using osmconvert (0.8.11) and by specifying a polygon boundary:
osmconvert large.osm.pbf -B=boundaries.poly -o=small.pbf

However, I'm getting the following error:
Segmentation fault (core dumped)

The poly file I'm using is retrieved from: https://download.geofabrik.de/europe/austria.poly
Unfortunately, nothing more I can share. No further details on the error itself. Is there a way to get smarter on what exactly could be going wrong? For example extra debug logging?


texts: In ArcMap, when using Add Surface Information, why are all avg_slope values positive?

texts:
I've got a polyline shapefile and a LIDAR elevation dataset.  I've successfully used Add Surface Information (from 3D Analyst) to calculate the average slope (avg_slope) of each line segment.  The problem is that all values are positive, no matter if the line is heading downhill or uphill.  I've checked the geometry of the line segments to ensure that the starting points of the segments I believe are heading downhill are indeed at a higher elevation than the end points...and they are.  When building a contour graph using the spatial analyst, it tells me the slope is downhill.

I'm working with a large set of lines, and I don't want to manually check each segment to see if it should have a positive (uphill) or negative (downhill) slope value.  So,

why doesn't this tool assign negative values for downhill lines, and
what is an automatable process that will compute downhill segments as having negative slope values?

An alternative way of stating the question is: if a polyline had a segment running from point A (elevation 200 ft) to point B (elevation 100 ft) and another segment running back from point B to point A, the Add Surface Information tool would return the exact same value for avg_slope, when one should be positive (B to A, going uphill) and one should be negative (A to B, going downhill).

How do I achieve this result?


texts: QGIS re-write loop efficiently, assign number to every row of a column,

texts:
I have a column titled CompanyNam and I need to extract the ownership information. For example, if CEZ is the first company then, I assign to it number 1, then Sokolovska is the next number, etc. If CEZ appears again in some other row, it will get number 1. It is important to note that if I have NULL in the column I assign a different number to every NULL row entry. I need the numbers corresponding to the CompanyNam output. I have the following code:

EmptyArray = []
d = {}
newlist = []
for gFeat in GeneratorLayer.getFeatures():
    Owner = gFeat.attributes()[gProvider.fieldNameIndex('CompanyNam')].toString()
    A = ([str(i) for i in Owner]) #convert from PyQt4.QtCore.QString to normal string
    B = ''.join(A)
    EmptyArray.append(B)
    for m, n in enumerate(EmptyArray):
        if n not in d:
            d[n] = [m+1]
        newlist.append({n: d[n]})
        if n == '':
            d[n] = [m+1]
        newlist.append({n: d[n]}) #Every NULL gets a new number
    for names in newlist:
        for o, p in names.iteritems():
            if o == '':
                a2 = str('{},NULL'.format(p))
            elif o != '':
                a2 = str('{},{}'.format(p,o)

I then use a2 in further steps. The code runs fine for a column with 60-100 rows but for larger columns, the computational times are really high. Can you suggest any way in which I can re-write this code, keeping the logic?
The output looks like this:

[1],CEZ
[1],CEZ
[1],CEZ
[1],CEZ
[1],CEZ
[1],CEZ
[1],CEZ
[1],CEZ
[9],Sokolovska
[10],International
[11],ENERGOTRANS,
[12],Alpiq
[13],Mittal Steel
[14],United
[1],CEZ
[1],CEZ
[17],Dalkia.....

It would be even nicer to have the numbering [1 ], [2], [3] and not [1 ], [9], [10], but I haven't figure out how to do that.


texts:
I am working with QGIS and I want to connect with Excel. I have exported attribute tables to Excel. Now I want to connect Excel with QGIS in the way that if I change a line dimension, it will be reflected in both the attribute table and in the Excel file.

Is there a way?


texts: I'm trying to plot spend data on zipcodes but only have US domestic zipcodes. Where can I find a shapefile/feature class for international postal codes? Point data is preferred but polygons will work.


texts: How to convert GeoTiff to PNG + World file in QGIS

texts:
I have quite a big area with DTM as a GeoTIFF. I want to save it to PNG with world file, but how to do it just in QGIS?

I know I can export view to PNG with world file, but it looses quality (Project/Export) .
I know how to transform format from TIFF to PNG but there I can't find any option for world file creation...
While i try to specify "PGW + yes" as extra option it says it is not supported...

What did I miss?


texts:
This is my first experience of using GIS systems.
Firstly, to give you an overview of my project, I am trying to analyze the obesity rates within the boroughs of London. I have gotten to the point now whereby I am trying to buffer the attribute of fast food restaurants whilst intersecting it with schools.
Unfortunately my unit of measurement for the buffer tool on both schools and fast food restaurants is in degrees. I have tried changing the CRS of both individual layers to a meters based CRS, but the layer just disappears when I do this and won't come back until I revert it to its previous CRS.
Why does that happen?


texts: Filter Image Collection by Multiple WRS Path/Rows using Google Earth Engine JavaScript API

texts:
How can I filter an image collection in Google Earth Engine by multiple WRS paths/rows?

This is an example of a filter using one path.row:

var sr2011 = ee.ImageCollection('LANDSAT/LT5_SR')
  .filterDate('2011-06-20', '2011-08-31')
  .filter(ee.Filter.eq('WRS_PATH', 59))
  .filter(ee.Filter.eq('WRS_ROW', 18));


texts: ol.interaction.Draw deleteCondition not working with Point feature

texts:
I'm not able to delete a point no matter which condition I write in deleteCondition.
It works with LineString and Polygon but not with Point features.

My version of OpenLayers is 3.19.1

Here is an excerpt of the code I'm using (it's TypeScript in Angular2 btw):

private _drawnFeatures: ol.Collection<any> = new ol.Collection<any>();

this._drawingSource = new ol.source.Vector({ features: this._drawnFeatures, wrapX: false });

this._drawingLayer = new ol.layer.Vector({
    source: this._drawingSource
});

this._drawInteraction = new ol.interaction.Draw({
    features: this._drawnFeatures,
    type: type,
    freehandCondition: ol.events.condition.altKeyOnly
});

var modify = new ol.interaction.Modify({
    features: this._drawnFeatures,
    deleteCondition: evt => { return ol.events.condition.shiftKeyOnly(evt) && ol.events.condition.singleClick(evt) }
});


texts: Is georeference the same as geotransform?

texts: Does georeference mean the same as geotransform?  If not then what is the difference?


texts:
How do I find the lat/lon of each pixel inside a landsat scene.
Either the equations for converting i,j to lat/lon or
just a bunch of files, one for each scene.
Background: I have downloaded a bunch of scenes from Greenland and now I want to
compare them to some MODIS data.  MODIS provides files with lat/lon for every pixel.


texts: Seeking arcpy.SetParameterAsText equivalent in ModelBuilder?

texts:
I know it sounds weird, but I would need a tool made with ModelBuilder that output two parameters:

a layer
a message (a string)

The constraint derives from the need to share it as a Geoprocessing service which a web application can use, but this application requires the two aforementioned outputs.

I know I can do it using Python, but I am explicitly asking if this is doable from ModelBuilder as well.

Basically, what I would like to translate from arcpy to ModelBuilder tool is:

arcpy.SetParameterAsText(0, my_layer)
arcpy.SetParameterAsText(1, my_message)

What I have now is a super simple MB tool like the one below:

However, I am not aware of any method to set the "String" variable as my output.

Using @Hornbydd’s suggestion, I finally managed to make it work using the Calculate Value tool in ModelBuilder.

The way I did it is shown in the picture below:


texts: Multiplying points to another feature table using ArcPy?

texts:
I have a point feature class like this:

And another table looks like this:

How can I add point's xy coordinate to that table based on their No. and automatic multiply as their No. multiplies using ArcPy?
In the end, it should just looks like this:


texts:
When wanting to use the different modules in GRASS 7.2.1 I get this error. The version of qgis is 2.14.17 on 64 bit windows.

I already checked that grass can find the right provider folder path.

Who knows how to get grass working?


texts: QgsSingleBandPseudoColorRenderer not automatically applied after added

texts:
Summary

I am trying to apply a Single Band Pseudo Color Renderer (SBPCR) to a raster layer, importing the shader from an XML. I've hit a dead end trying to get my SBPCR to apply.

QGIS Version - 3.4.5
Mac OS X v10.14.4
Python 3.6

Problem

Unlike my hillshade raster which, upon being added to the renderer seems to refresh automatically, this does not happen with my SBPCR renderer.

Code Example

def add_color_gradient_raster(color_gradient_rlayer):
    # Import Color Ramp
    my_qfile = QFile("/Hydro1/top_15lev.xml.xml")
    my_qdom = QDomDocument()
    my_qdom.setContent(my_qfile, False)
    my_raster_shader = QgsRasterShader()
    my_raster_shader.readXml(my_qdom.documentElement())
    # Single Band Pseudo Color Render Created
    sb_pseudo_color_renderer = QgsSingleBandPseudoColorRenderer(
    color_gradient_rlayer.dataProvider(),
    1,
    my_raster_shader)
    sb_pseudo_color_renderer.setOpacity(0.6)
    color_gradient_rlayer.setRenderer(sb_pseudo_color_renderer)
    # Add layer into project
    QgsProject.instance().addMapLayer(color_gradient_rlayer, False)
    layerTree = iface.layerTreeCanvasBridge().rootGroup()
    layerTree.insertChildNode(3, QgsLayerTreeLayer(color_gradient_rlayer))

Output

As seen above I have reason to believe that the problem is not in the creation of the renderer.

The layer is successfully created and the renderer's opacity is set.
Upon editing the layer's properties all options seem to be set
correctly (including the imported colours) and I only need to hit
apply/ok.

How do I "apply" these settings?

I am not sure what I am missing.


texts: Advantages of R-trees in comparison to geohashes

texts:
Geohashes are being widely used in products like: Lucene, mongodb, etc and have become one of the most important technology of present day.

Have Geohashes replaced the good old R-trees or do R-trees have any advantages as compared to Geohashes?


texts: Calculating median value of polygon from points in ArcGIS Desktop?

texts:
I have a point layer coming from LIDAR data containing x,y,z and intensity values. On the other hand I have a layer containing several polygons. All the points are inside the polygons.

I need to calculate the median Z value of each polygon using the points. I could do a spatial join, but it offers only the possibility of calculating mean values unfortunately, not median.

How do I proceed?


texts: Can you draw an actual line instead of a cutline with GDAL?

texts:
I'm fairly new to GDAL and I have some GeoTIFF files that I would like to highlight some bodies of water on.  I have gotten pretty far, but I have one problem.  The only way I am able to highlight a body of water is with gdal.warp and cutlineDSName (I'm using osgeo gdal in Python), which is great, but it cuts out the rest of the image completely.  What I would like to do is just draw a line, rather than cut the surrounding area out.

Can anyone help me out?

I'm hoping it's a flag that I am missing somewhere.

code below.

tifFile = 'C:\\merged.tif'
shapeFile = 'C:\shapefile.shp'
outFile = 'C:\\merged_and_cropped.tif'

input_raster = tifFile
output_raster = outFile
input_kml = shapeFile

ds = gdal.Warp(output_raster,
              input_raster,
              format = 'GTiff',
              cutlineDSName = shapeFile,
              cropToCutline=True,
              dstNodata = 0)
ds = None


texts: Selection and Popups on Overlapping/Coincident Vector Layer Polygons in OpenLayers

texts:
I have a GeoJSON file of polygons that I am bringing into my OpenLayers session.  Some of the polygons occupy the same space on the map.  I am using a select control and popup control like:-

selectControl = new OpenLayers.Control.SelectFeature(geojson_layer,
{
    onSelect: onFeatureSelect, onUnselect: onFeatureUnselect,
    clickout: true, toggle: true,
    multiple: false, hover: false,
    toggleKey: "ctrlKey", // ctrl key removes from selection
    multipleKey: "shiftKey" // shift key adds to selection
});
map.addControl(selectControl);
selectControl.activate();
geojson_layer.events.on(
{
    "featureselected": function(e)
    {
        showStatus("selected feature "+e.feature.attributes.APPREF+" on  Applications");
    },
    "featureunselected": function(e)
        {
        showStatus("");
        }
    });
        function onFeatureSelect(feature)
        {
            selectedFeature = feature;
            popup = new OpenLayers.Popup.AnchoredBubble
            (
                "chicken",
                feature.geometry.getBounds().getCenterLonLat(),
                null,
                "
App Reference: " + feature.attributes.APPREF +"
Proposal:
" + feature.attributes.PROPOSAL+ '
Further Details'+"
",
                null,
                true,
                onPopupClose
            );
            feature.popup = popup;
            popup.autoSize = true;
            popup.minSize = new OpenLayers.Size(150,50);
            popup.maxSize = new OpenLayers.Size(250,150);
            //popup.panMapIfOutOfView = true;
            popup.closeOnMove = true;
            map.addPopup(popup);
        }
        function onFeatureUnselect(feature)
        {
            map.removePopup(feature.popup);
            feature.popup.destroy();
            feature.popup = null;
        }

Where more than 1 polygon occupies the same space - can anyone advise how I allow for drilling down and/or making the other features at the location available for a popup (and therefore the link out for further information).

I suspect the answer lies somewhere in Z-indexing, but can't wrap my head around how I would integrate that from the example thats provided.

Presumably there is a way of handling this type of scenario - if someone could enlighten me that would be very much appreciated

Thanks


texts: Creating geometry from lat/lon in table using PostGIS?

texts: I have a table with over 800,000 records including a lat and long column. The data source did not create geometries, so I was looking around how to do this. I know there's ST_MakePoint(x,y,SRID) but all the documentation and other threads show how to create single points. How do I create the geometry type using these columns?


texts: How do I get vector feature data from an OpenStreetMap layer?

texts:
I imported the OSM tile layer and can see it visually. But how do I access the actual data of the layer? It's vector data, made of points, lines, and polygons, right? How can I access the feature data for some part of the map, like as json?

Ex: Get the data for the polylines/polygons making up the streets/blocks of a town.


texts: Disable warning that scratch layer isn't saved?

texts:
I'm upgrading a QGIS plugin that I've written for farmers (GeoDataFarm) and I wounder if it is possible to disable the warning of temporary scratch layers on close project that was added in QGIS 3.2?

I use the scratch layer for the user to draw the field on the map, then I use the polygon in wkt format to store the data on a Postgres database.
A more elegant solution might be to remove my memory layer and add it again. Hmm, any thoughts?

Here is some of my code (only some tiny bits to make sense of my question).

def add_field(self):
    self.field = QgsVectorLayer("Polygon?crs=epsg:4326", name, "memory")
    self.field.startEditing()
    self.iface.actionAddFeature().trigger()
    QgsProject.instance().addMapLayer(self.field)

def save_field(self):
    self.iface.actionSaveActiveLayerEdits().trigger()
    self.iface.actionToggleEditing().trigger()
    feature = self.field.getFeature(1)
    polygon = feature.geometry().asWkt()
    self.db.execute_sql('insert into fields (geom) values (st_geomfromtext('{wkt}', 4326)'.format(wkt=polygon))


texts: How to Disable the add-in button when ArcMap is started?

texts:
Question 1:How to Disable the add-in button when the ArcMap started,then I use the SelectByRectangletool  select features,then the button enabled! I am a fresh man about the add-in.

Question 2:I made a menu as show ,As the figure shows, how to avoid the part of in red rectangle,because the menu do not have the child menu again,but it is show the arrow ,and part of the child menu.

I means the tool’s skin is gray，and it is not activate

I am sorry about my poor engilsh

here is the question2's config.esriaddinx file:

<ArcMap>
  <Commands>
    <Button id="Yangtze_University_GRP_Mult" class="Mult" message="new" caption="test" tip="tip" category="Add-In Controls" image="Images\1.png" />
    <MultiItem id="ESRI_SelectionSample_MultiItem" caption="MultiItem AddIn" class="ZoomToLayerMultiItem" >
    </MultiItem>
  </Commands>
  <Menus>
    <!--main menu-->
    <Menu id="Yangtze_University_GRP_My_Menu" caption="MainMenu" isRootMenu="true">
      <Items>
        <Menu refID="Yangtze_University_GRP_My_Menu_2" />
        <Menu refID="Yangtze_University_GRP_My_Menu_1" />
      </Items>
    </Menu>
    <!--submenu-->
    <Menu id="Yangtze_University_GRP_My_Menu_2" caption="New file" isRootMenu="false">
      <Items>
        <!--<Button refID="Yangtze_University_GRP_Mult" />-->
        <Menu refID="Yangtze_University_GRP_My_Menu_3" />
        <Menu refID="Yangtze_University_GRP_My_Menu_4" />
        <Menu refID="Yangtze_University_GRP_My_Menu_5" />
        <MultiItem refID="ESRI_SelectionSample_MultiItem" />
      </Items>
    </Menu>
    <!--third menu-->
    <Menu id="Yangtze_University_GRP_My_Menu_3" caption="test1" isRootMenu="false" isShortcutMenu="true" />
    <Menu id="Yangtze_University_GRP_My_Menu_4" caption="test2" isRootMenu="false" isShortcutMenu="true" />
    <Menu id="Yangtze_University_GRP_My_Menu_5" caption="test3" isRootMenu="false" isShortcutMenu="true" />
  </Menus>
  <Toolbars>
    <Toolbar id="Yangtze_University_GRP_My_Toolbar" caption="toolbar" showInitially="true">
      <Items>
        <Menu refID="Yangtze_University_GRP_My_Menu" />
      </Items>
    </Toolbar>
  </Toolbars>
</ArcMap>


texts:
I've searched a lot but didn't find an answer: When editing lines with QGIS, I often need to split a linestring exactly on one of his existing internal vertices, but I can't find a way.
I've tried the tool "Split Features" thinking that it could snap on vertices, but I was wrong, and I've tried several ways of selecting the vertex with no luck... Any solution?


texts: Calculate Row if Row Value Equals Value in List

texts:
I have a folder of photos named with a building ID that corresponds to a building ID in a building footprints feature class. Not all rows in the footprints feature class have a corresponding photo. So, I want to use an update cursor to calculate a photofile field in the footprints feature class if there is a photo in the folder.

I've gotten this to work with both arcpy.UpdateCursor [getValue(), setValue()] and arcpy.da.UpdateCursor, but I don't think it's the most efficient method due to the time it takes to run (new cursor for each value in the list?). Here is the working method using da, but again, it takes forever to run.

import arcpy

photos = arcpy.GetParameterAsText(0) #photos folder
fc = arcpy.GetParameterAsText(1) #feature class

arcpy.env.workspace = photos
photoList = arcpy.ListRasters()

fields = ('BLDG_ID','PHOTOFILE')

with arcpy.da.UpdateCursor(fc,fields) as uc:
    for row in uc:
        for photo in photoList:
            if (row[0] == photo[:-4]):
                row[1] = photo
            else:
                pass
            uc.updateRow(row)
del row, uc

That's probably not the most Pythonic way. I want to use something like:

with arcpy.da.UpdateCursor(fc,fields) as uc:
    for row in uc:
        if row[0] in photoList:
            row[1] = value in list
        uc.updateRow(row)
del row, uc

The problem is that the values in the photoList have ".jpg" after the BLDG_ID, so it is not equal without removing ".jpg" from each value in the list or concatenating ".jpg" to the BLDG_ID value in the row (row[0]+".jpg" ... I don't think that works). I'm also not sure how to update the row in the field as the value in the list without iterating through the list, and I need to retain the ".jpg" in the PHOTOFILE field as well.

Example:

Photo in folder = 12345.jpg

BLDG_ID in fc = 12345

PHOTOFILE in fc = 12345.jpg


texts: Selecting polygons not covered by any other polygon (or BBox) in table in PostGIS?

texts:
In a table where many of the geometries might be completely covered by other ones, is there a query that selects only geometries that are not covered by any other geometry in the same table (excluding the geometry itself of course)?

Each geometry has properties that should be returned too, so no new geometry should be created.


texts: What is the most efficient way to do a spatial join of all raster points within a polygon layer area?

texts:
I have a raster GeoTIFF of a digital elevation model and a polygon shapefile with spatial boundaries of the raster. Here is a snapshot of raster and polygon layer:

The shapefile has one attribute named "plot #".  My ultimate goal is to extract individual pixel values in the raster that fall under the polygons and join them onto the new shapefile that has plot# field. Currently, my workflow does this:
In Qgis, export raster to point file (huge file) >> perform spatial join over the plot boundaries using custom programs in R. The final output table should look like (bonus if it is a shapefile!):

Plot#  height

1 60

1 30

1 40

1 .

1 .

2 20

2 30

. .

Now the above workflow has a few issues. First, to me, the raster-to-point file export seems unnecessary and a slow process. Second, the spatial join is painfully slow for 1800 polygons on a dataframe. I wonder if:
1. there is any way to bypass the raster-to-point conversion step and directly apply the spatial join on raster pixels similar to what 'zonal statistics' plugin in this does but for all individual raster points within the polygon. I find zonal statistics to be very efficient and fast but it outputs only summary stats, not individual points.
2. If there is no other alternative, is it possible to do the same steps i.e. raster-to-point >> spatial join all in R instead of going back and forth in this and R?

How should I tackle the issue?


texts:
On occasionally when I open a raster in ERDAS 2013, I get white or yellow circle icons on my ERDAS Contents and tried to find information on ERDAS.

I want to know what does the icon mean ? Thanks appreicate it.

Her is an attachment of what it looks like


texts: Waiting for a python script tool to complete while using pythonaddin button

texts:
I have created a pythonaddin toolbar, with a button that executes a python script tool.  This tool is in an ArcToolbox which has been stored in the Install folder of the addin.  It executes and runs perfectly.  However, the rest of the script in the pythonaddin button class does not wait for the tool to finish.  I searched on the ESRI help forums, and there was a suggestion that I use a while loop to wait for the result of the tool.  However, I do not believe that a custom python script tool returns a result object like arcpy gp tools.

My toolbar (toolbar_addin.py) code

parameterList = []
class button1(object):
    def__init__(self):
        self.enabled=True
        self.checked=False
    def onClick(self):
        ### Get path to external toolbox
        toolbox = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'Toolbox.tbx')
        ### execute external script tool
        pythonaddins.GPToolDialog(toolbox, 'tool')

        ### wait for result from gp tool
        while len(parameterList) == 0:
            time.sleep(0.2)

        #do something with items from parameterList

My tool (tool.py) code

import toolbar_addin

input1 = arcpy.GetParameterAsText(0)
input2 = arcpy.GetParameterAsText(1)

#do something to create output1 and output2...

toolbar_addin.parameterList = [output1, output2]

Now I get the following error:

Traceback (most recent call last):
File "C:\path\to\addin\toolbar_addin.py", line 229, in onClick
    time.sleep(0.2)
TypeError: GPToolDialog() takes at most 1 argument (2 given)

Any suggestions on how I can force my code to wait on the gp tool?

EDIT: I've modified the fuctioning of my button and script tool to resemble this post.  The tool passes parameters to a global list variable "parameterList" in the toolbar.


texts: Translating .bil files in batch mode with QGIS?

texts:
I have some .bil files that I would like to convert into ASCII Gridded XYZ format files.

It seems to work well with the Raster>Translate tool of QGIS when I process one file at a time. But I would like to do it by batch. So I checked the batch box but got this error message:

Processing of the following files ended with error:

    C:/Users/Fall Line Tower/Desktop/Bil To Be Processed/Only Bill/PRISM_tmax_stable_4kmD1_19810101_bil.bil
    ERROR 4: `C:/Users/Fall Line Tower/Desktop/Bil To Be Processed/Onli Bill/PRISM_tmax_stable_4kmD1_19810101_bil.bil' not recognised as a supported file format.
    GDALOpen failed - 4

I have the same error message for each .bil file within my source directory.

Any idea where the problem comes from?

Does the batch mode take only one type of format?

Not the .bil?

I am using QGIS 2.2 on Windows 7.


texts: Calculate time since event in a raster time series in R

texts:
I have a time series of binomial data and I want to know the number of time steps since the last event (1). I am just not sure how to do this in R.
Here is an example data set. For each cell in the raster, I want to know the number of timesteps since a 1 was recorded.
r <- raster(nrow=100, ncol=100)
r[] = rbinom(10000,1,.1)

for (i in 1:9){
  tmp = raster(nrow=100, ncol=100)
  tmp[] = rbinom(10000,1,.1)
  r = stack(r,tmp)
}

r = brick(r)

So for each cell of the raster brick, you would have a series (e.g. c(0,0,0,0,1,0,0,0,1,0)). In this case the time to the event would by 4. I just can't figure out how to do this in a function that can be applied to the raster brick.


texts: Reading and parsing .grc/.grd files in Python

texts:
I am trying to read .grc/.grd file and convert the data into dataframe and write to Hive. One of our colleague did this in "R" and I am trying to convert this into Python.
I am using rasterio library to read the file but did not find any API like rasterToPoints to convert into dataframe
Please anyone who can help me to convert this into python


texts:
When importing a CSV file into Google Earth Pro, I am asked if I want to apply a style template. I do. I'd like to be able to load a previously defined template rather than have to set the style options manually each time I load a file. After setting all of the options in the "Style Template Options" dialog, it appears to save two styling-related files: template_name.kdx and template_name.kst.

My template is fairly involved with a lot of subfolders, so manually setting it up in the dialog box every time is tedious. It looks like the basic settings are recorded in the .kst file. How do I apply a saved template to my imported records?

I would also like to set icon size (or is it scale) in the template. Currently, I can only control the icon image/href used. I would also like to set the label in the template. I have been unable to find docs/specs for the .kst format. Are there additional configuration options that can be manually edited in the .kst file that are not configurable from the dialog box?


texts: How to assign colors to specific raster values using R

texts:
I am using the levelplot function of the awesome rasterVis package to create a map whose values diverge around zero. I am plotting a raster file using a red to blue palette, but I am trying to assign grey to zero values in the map.

Specifically, I want to reproduce the colors of this figure:

Notice that where a red to blue scale is used, but zero values have been colored with grey.

Currently, my maps looks like this:

And this is the code to reproduce it (file available at https://www.dropbox.com/s/cypfdu1eaz2fuok/r.annual.tif?dl=0):

# Load required packages
library(rasterVis)

# open file
r.annual <- raster("Downloads/r.annual.tif")

# Set color palette
myTheme=rasterTheme(region=brewer.pal('RdBu', n=11))

# Plot
levelplot(annual.mask, par.settings=myTheme, margin=F)

How can I assign grey to all zero values in my map above?


texts:
I have a model that generates a list of materials for a proposed project. It has no parameters which is what I want, and it runs almost perfectly. For some reason if I open it to EDIT and run the process it generates my .XLS file perfectly, but if I just run it from Toolbox it'll generate the file and quickly delete it.

Does anyone know why this happens?


texts: Merging two spatial objects in R: spCbind gives error “row names not identical”

texts:
## Importing map of Brazil
map <- readOGR ("UFEBRASIL.shp")
class (map)
plot (map, axes = F)
title ("")

## importing file with case data
data <- read.csv2 ("cases.csv")

# Recognizing the names of states as a character
map @ data $ NM_ESTADO <- as.character (map @ data $ NM_ESTADO)
map @ data $ NM_ESTADO

# Ordering the data correctly to match the order of the neighborhoods in the shape and attribute table.
head (map @ data)
head
ind <- match (map @ data $ ID, data $ ID)
ind
# Correcting accents in shape
proj4string (map) <- CRS ("+ proj = longlat + datum = WGS84 + no_defs")
Encoding (map $ NM_ESTADO) <- "UTF-8"
map $ CD_GEOCODU [is.na (map $ CD_GEOCODU)] <- 0
map1 <- spTransform (map, CRS ("+ proj = longlat + ellps = GRS80"))
summary (map1)
map1 $ NM_ESTADO

# Concatenating the map with the database
data <- data [ind,]
row.names (data) <- map1 $ NM_ESTADO
map1 <- spCbind (map1, data)

In this last line of command of the following error:

  Error in spCbind(mapa1, dados) : row names not identical


texts:
I know people have already asked this question before but the last time was about 4 years ago so I was wondering if something changed since.
So, does anyone know how to perform kriging in qgis? and not via SAGA but actually in the program itself.


texts:
We are having a client which is looking to develop a service provider directory , which will be provided free of cost to end users nothing commercial for end users they can use it with out any charges.

End users will be able to search providers on the directory and result of search will be show on google maps also end users will be able to  search directory pointing different locations on map so directory will show results for that particular location itself.

Everything for End user is going to be free. But the service providers are going to be charged for their subscriptions and that will be minimal plus some part of this charges will go do different charities.

So question is that do they have to take Google Enterprise version? Or Google free version will do?

Thanks,
Naresh.


texts: Leaflet Routing Machine: how to export route coordinates with time associate

texts:
I'm using Leaflet Routing Machine with Graphhopper as router. I just want to get routes coordinates and time associate to them. I don't know if it's possible. I only have summary.totalTime and instructions.times.

But what i want is time according to latLng

        // Route informations
        track.on('routeselected', function(e) {
            var route = e.route;
            var rCoordinates = [],
                i,
                latLng;

            for (i = 0; i < route.coordinates.length; i++) {
                    latLng = L.latLng(route.coordinates[i]);
                    rCoordinates.push([latLng.lng, latLng.lat]);
            }

        console.log(route.summary.totalTime);
        console.log(route.coordinates.length);
        console.log(route);


texts: K-means clustering in QGIS with statistically significant difference

texts:
Anyone knows that how to do K-means clustering with significant difference?
What I means is as the image shows, the area is divided into two zones using k-means clustering, and the significant difference is also obtained.
I found tools that can do k-means clustering, but I have no idea how the p-value is calculated.


texts: Generating DEM from contour Shapefile using QGIS?

texts:
I want to create a DEM from depth contours (shapefile) in QGIS. I've been trying to work out how to do this based on this previous question, notably from this link from that page, however I'm stuck, and was hoping someone could help.

My progress/process so far:

open new map in  QGIS, import the contour shapefile
set GRASS extent, set new GRASS mapfile, etc
try to 'import' the shapefile into grass by converting it into a grass file, based on a hunch from the wording in this page (search in page: 'use of qgis and grass'): file management, import into grass, import vector into grass, import loaded vector (v.in.ogr.qgis), then view output. This creates an identical looking layer.
try to use r.surf.contour on this - doesn't allow layer to be selected.
try to use r.surf.contour2 on this - either crashes QGIS or produces a full block of colour.

Any thoughts/tips on what I should be doing differently?

Interpolation failing, tried everything in GRASS & in GRASS plugin, failed, could problem be due to source shapefile?

Something I learnt from elsewhere was 'geometry tools > check geometry validity' which I ran on the contour shapefile and it found loads of errors.

I'm wondering if this could be the reason that all my efforts to create a DEM are failing?

Would anyone do me the huge favour of having a look at my shapefile to see if I'm just being useless, or if it's a problem that a relative novice would be excused for not anticipating?


texts: ExecuteError: ERROR 999999 with "PointsToLine" function with a Point Feature class created by "MakeRouteEventLayer_lr"

texts:
I have an issue with the "PointsToLine" function in Python arcpy and in the tool dialog, the problem comes everytimes I use a POINT Feature Class (FC) created from the "Make Route Event Layer" function (Linear Referencing).

All the properties seem to be the same as a regular (working) POINT class, but the error keeps prompting with such created FC. Is there something I am missing? Is this a known issue? Is there a workaround for this?

Here is my simplified code to generate the Point FC used for the "PointToLine" function:

import arcpy

table_points_TEST = "C:\\ArcGIS\\Test.gdb\\TEST_TABLE_TO_LINE"
axe_repere = "C:\\ArcGIS\\Test.gdb\\AXE_REPERE"
table_points_TEST_TEMP = "C:\\ArcGIS\\Test.gdb\\TEST_TABLE_TO_LINE_TEMP"
FC_point = "C:\\ArcGIS\\Test.gdb\\FC_POINT"
FC_line = "C:\\ArcGIS\\Test.gdb\\FC_LINE"

arcpy.MakeRouteEventLayer_lr(axe_repere , "AXIS_NAME", table_points_TEST, "AXE POINT M_POINT", table_points_TEST_TEMP)
arcpy.CopyFeatures_management(table_points_TEST_TEMP, FC_point)

arcpy.PointsToLine_management(FC_point,FC_line)

If I use the Point layer generated by the "CopyFeatures" in the "PointToLine" tool dialog, the same  error occurs.


texts: How to find least-cost paths including tunnels (where slopes are too steep)?

texts: I'm currently working on a road path finding project for a motorway in ArcGIS 10. One criterion is obviously the slope angle. What I would like to do is to create a new cost surface which finds possible crossings where high mountains are --> Where would be a good path for a tunnel instead of going around the mountain?  Any ideas how I could do that?


texts: QFields often don't start-> error "QField angehalten"

texts:
I run QField on Android Samsung S7 with one small project. Sometimes it runs ok but actually I can't open the program. When I try it several times I get the message "QField bleibt ständig stehen/QField stops permanently" after QField "QField angehalten/QField stopped". I used the version from the GoogleAppstore 1.6.4 and also after deinstallation I tried the APK version 1.7.3 (qfield-v1.7.3-armv7.apk)? No improvement.
Does anyone have the same error pattern?


texts: Watershed tool producing same output as basin tool for all pour points

texts:
I am using the Hydrology Tools in the Spatial Analyst Toolbox (ArcGIS 10.3)

I am trying to delineate watersheds in a region in order to do spatial interpolation by top-kriging at some ungaged sites.  I have followed all the intermediary  steps for the delineation (Filling, Flow Direction, Flow Accumulation, Pour Point Snapping) successfully. However, when I get to the Watershed tool, regardless of where I place the pour points, the output I get is the basin (generated with the Basin tool) in which the pour point is located. For example in the image attached, the output from the basin tool is shown on the left. When I run the watershed tool with the snapped pour point in red, the delineated watershed (in green, on the right) is just the basin in which the pour point is located instead of the upstream contributing catchment for that particular pour point.

Does anyone know what I might be doing wrong?


texts: I am interested in knowing about current trends in application of GIS for medical geography other than disease/epidemic mapping, heath atlas, availability of health infrastructure (hospitals, blood banks, primary health centers etc).