texts: Vector layers without geometry in PyQGIS?

texts:
I have set up a couple of automated layer additions via the toolbar.
After adding geometry in each respective layer (point, lineString, polygon) and saving changes to the layers, the added forms are non-present in the layers upon restart and reopening of files.

Additionally, to actually get 'run3' to pass and the layer to be added, the button has to be pressed twice to proceed through my if statements.  Which I don't clearly understand.

Please take error and non-Pythonic ways in the following code as a product of my novice skillset.

    # run add base
def run3(self):
    if FindDirectory._destination is 0:
        projectDirectory = bool(QgsProject.instance().readPath(''))
        if projectDirectory:
            projectDirectory = QgsProject.instance().readPath("./")
            FindDirectory._destination = str(projectDirectory)
        else:
            iface.messageBar().pushMessage("Error", "Please select a central folder location to save your base layer to!", level=QgsMessageBar.CRITICAL, duration=3)
    elif FindDirectory._destination is not 0:
        if GetMaps._TOPOcrs is 0:
            canvas = bool(self.iface.mapCanvas())
            if canvas:
                canvas = self.iface.mapCanvas()
                mapRenderer = canvas.mapRenderer()
                srs = mapRenderer.destinationCrs()
                GetMaps._TOPOcrs = str(srs.authid())
            else:
                iface.messageBar().pushMessage("Error", "Please select a topographic map for your simulation!", level=QgsMessageBar.CRITICAL, duration=3)
        elif GetMaps._TOPOcrs is not 0:
            crs = str(GetMaps._TOPOcrs)
            Base =  QgsVectorLayer('Polygon?crs=' + GetMaps._TOPOcrs + '&field=BOTTOM:string(6)&field=TOP:string(6)&field=K:string(8)&field=Porosity:string(6)&field=Recharge:string(8)&field=IHead:string(5)&field=IBound:string(3)&index=yes', "Base", "memory")
            QgsMapLayerRegistry.instance().addMapLayers([Base])
            pr = Base.dataProvider()
            BaseFeat = FindDirectory._destination+"\Base.shp"
            writer = QgsVectorFileWriter(BaseFeat, "CP1250", pr.fields(), QGis.WKBPolygon, pr.crs(), "ESRI Shapefile" )
            iface.setActiveLayer(Base)


texts: PostgreSQL PL/Python Function for "IF" condition

texts:
I am using OSM data to run pgr_dijkstra() for weighted-path calculation. But after importing OSM data using osm2pgsql I populated some additional columns into my ways table. So now it has gid serial, the_geom geometry, source integer, target integer, source_height double precision, target_height double precision plus some other columns. The source_height and target_height columns I have populated using the DEM of the region.

Simply running the pgr_dijkstra() on ways will give seq, id1 (node), id2 (edge), cost columns. What I want is to also determine the source_height or target_height for a particular edge if the node is equal to the source or target. I was trying to write a procedural function using PL/Python (although I don't know which procedural language is better) but couldn't complete.

Following is the attempt to create a PL/Python function. The problem is I don't know how to append all the rows coming from the for loop into a variable to return it as a table.

CREATE OR REPLACE FUNCTION test1(source int, target int, tablename text)
RETURNS setof varchar # I am not sure if this return type is correct
AS
$$
import psycopg2
import osgeo.ogr
statement1 = ("select count(*) from pgr_dijkstra(' SELECT gid AS id, source::integer, target::integer, length_dem::double precision AS cost FROM %s', %s, %s, false, false) a LEFT JOIN %s b ON (a.id2 = b.gid)" % (tablename, source, target, tablename))
statement2 = ("select seq, id1 AS node, id2 AS edge, cost, b.the_geom, b.source, b.target, b.slope_st_pt, b.slope_end_pt from pgr_dijkstra(' SELECT gid AS id, source::integer, target::integer, length_dem::double precision AS cost FROM %s', %s, %s, false, false) a LEFT JOIN %s b ON (a.id2 = b.gid)" % (tablename, source, target, tablename))
count1 = plpy.execute(statement1)
run2 = plpy.execute(statement2)
len1 = count1[0]['count']
row1 = [];
for i in range(len1):
    if run2[i]['node']==run2[i]['source']:
        # Code to append node, edge, the_geom and slope_st_pt as slope to row1
elif run2[i]['node']==run2[i]['target']:
        # Code to append node, edge, the_geom and slope_end_pt as slope to row1
return row1; # Return row1 table just created using the for loop above
$$
LANGUAGE 'plpythonu' VOLATILE;


texts:
I try to create MapView without XAML but I can't. My code blocks in below

class ShowMap : ContentPage
{
    public ShowMap()
    {
        Title = "Show a Map";

        Map myMap = new Map(Basemap.CreateImagery());
        MapView MyMapView = new MapView();
        MyMapView.Map = myMap;

        //Content = new StackLayout
        //{
        //    VerticalOptions = LayoutOptions.CenterAndExpand,
        //    HorizontalOptions = LayoutOptions.CenterAndExpand,
        //    Children =
        //    {
        //        MyMapView
        //    }
        //};

    }
}

If I uncomment block in code. Shown an error.

cannot convert from 'Esri.ArcGISRuntime.UI.Controls.MapView' to 'Xamarin.Forms.View'

How can I do this?


texts:
I'm working on QGIS 2.18.13 version on my Mac. Now I want to change to 2.18.15 version. On Windows, it's really easy to make.

After uninstalling from the control panel, you go through the .qgis2 to remove it and it's completely gone. However, on Mac, I can't find this folder.

Could you please describe to me the method to follow to uninstall it?


texts: Where to download Print Widget for ArcGIS Flex Viewer 3.4?

texts:
I would like to download the Print Widget for ArcGIS Flex Viewer 3.4 which contains Scale Units like below:

Does anyone know where I can do that, please?


texts:
When a user clicks up an infowindow, it gets all the data that I set through the CartoDB InfoWindow UI.

Sometimes, this data is a bit too long to fit into a small infowindow, so I want to truncate that data to say a 100 characters.

How do I edit the data before it is shown to the user in an InfoWindow?
I'm creating the Vis(ualisation) from an existing visualisation on my CartoDB account using CartoDB.js.


texts: In the metadata for GIS shapefile, what is the syncdate?

texts: I'm trying to discern when the files were last updated, but the sync date appears to be more recent than the date created, so I just wanted to check if sync date was the date the file was updated?


texts: Can I hide the "Map created by {user}" element in a iframe visualization?

texts: I would like hide the "Map created by {user}" element that its show when you use iframe embed.


texts: Why is WPS execute process in Geoserver producing this error?

texts:
I am trying to compute raster zonal statistics using a TIFF and a Shapefile. I am using WPS process extension for GeoServer. I am trying to use the WPS request builder provided by my GeoServer installation under demos.

I am trying to execute the process to see the output but I am getting the message in the second screenshot. My geoTiff file is 459 MB. Am I doing something wrong?

The following are the screenshots:

.

The following is the xml output generated by clicking the "Generate xml from process I/O":

<?xml version="1.0" encoding="UTF-8"?><wps:Execute version="1.0.0" service="WPS" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.opengis.net/wps/1.0.0" xmlns:wfs="http://www.opengis.net/wfs" xmlns:wps="http://www.opengis.net/wps/1.0.0" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:wcs="http://www.opengis.net/wcs/1.1.1" xmlns:xlink="http://www.w3.org/1999/xlink" xsi:schemaLocation="http://www.opengis.net/wps/1.0.0 http://schemas.opengis.net/wps/1.0.0/wpsAll.xsd">
  <ows:Identifier>ras:RasterZonalStatistics</ows:Identifier>
  <wps:DataInputs>
    <wps:Input>
      <ows:Identifier>data</ows:Identifier>
      <wps:Reference mimeType="image/tiff" xlink:href="http://geoserver/wcs" method="POST">
        <wps:Body>
          <wcs:GetCoverage service="WCS" version="1.1.1">
            <ows:Identifier>soil_system:PI</ows:Identifier>
            <wcs:DomainSubset>
              <ows:BoundingBox crs="http://www.opengis.net/gml/srs/epsg.xml#4326">
                <ows:LowerCorner>-128.14691380352517 22.714877523496572</ows:LowerCorner>
                <ows:UpperCorner>-65.08558771879018 51.75931655065378</ows:UpperCorner>
              </ows:BoundingBox>
            </wcs:DomainSubset>
            <wcs:Output format="image/tiff"/>
          </wcs:GetCoverage>
        </wps:Body>
      </wps:Reference>
    </wps:Input>
    <wps:Input>
      <ows:Identifier>zones</ows:Identifier>
      <wps:Reference mimeType="text/xml" xlink:href="http://geoserver/wfs" method="POST">
        <wps:Body>
          <wfs:GetFeature service="WFS" version="1.0.0" outputFormat="GML2" xmlns:soil_system="soil_system">
            <wfs:Query typeName="soil_system:USA_States_WGS84"/>
          </wfs:GetFeature>
        </wps:Body>
      </wps:Reference>
    </wps:Input>
    <wps:Input>
      <ows:Identifier>classification</ows:Identifier>
      <wps:Reference mimeType="image/tiff" xlink:href="http://geoserver/wcs" method="POST">
        <wps:Body/>
      </wps:Reference>
    </wps:Input>
  </wps:DataInputs>
  <wps:ResponseForm>
    <wps:RawDataOutput mimeType="text/xml; subtype=wfs-collection/1.0">
      <ows:Identifier>statistics</ows:Identifier>
    </wps:RawDataOutput>
  </wps:ResponseForm>
</wps:Execute>


texts: PostGis Dump file to SpatialLite without going through PostgreSQL

texts: I received a large (2Gb) postgis dump file that I'd like to work with. But I don't want to deal with a PostgreSQL installation so I though about using SpatialLite. Is it possible to restore a pg_dump file straight into Sqlite?


texts: Add points to shapefile and display over map interactively using Python Add-Ins

texts:
I want to add points to a shapefile interactively using pythonaddins and display them over map. Points are added to shapefile but not to display. Please suggest how to fix this discrepancy.

class AddPoints(object):
   """Implementation for points_addin.tool1 (Tool)"""
   def __init__(self):
       self.enabled = False
       self.shape = 3

   def onMouseDownMap(self, x, y, button, shift):
       xy = (x,y)
       cursor = arcpy.da.InsertCursor(r'C:/temp/points.shp', ("SHAPE@XY"))
       cursor.insertRow([xy])

       arcpy.RefreshTOC()
       arcpy.RefreshActiveView()


texts: How to find the distance of a point from the begning of a StringLine with shapely?

texts:
I have a route (represented by StringLine) and a location (Point) and I need to find the  progress of the location in the route. The progress is the distance in meters from the begning of the route. The location will not necessarily be excactly on the line so I take a little buffer (lats say 30m).

This is what I got so far though after testing a bit I think this is not the correct solution

if point.buffer(float(radius_in_meters)/100000).intersects(route):
   distance = route.project(point) * 100000

I'm new to this stuff and i really hope i wont have to implement it my self.

P.S I have no sentiments for shapely so if there is a different recommanded package for python it's also fine.


texts:
I am looking to run MST on some vectors with geometries and cost. I have not been able to find a one stop python/arcpy method to do so. The closest I've gotten is by using the code in the comment by dkweins at:

https://community.esri.com/blogs/dan_patterson/2017/01/31/spanning-trees

This code does not have a place to assign weight to the graph.

Is there a package MST tool in Arc where I can just input a set of points and receive the paths out in consideration of cost.
If no (1), what can I do to the linked code to allow this functionality?

I have tried adding a cost field to the search cursor but I can't figure out how to append this info to the dictionary Q. If I could, the edit to addToTree should be pretty straight forward.

ORIGINAL:

def addToTree(F,Q,sr):
    min_dist = float("inf") # set to infinity
    for fk,fv in F.iteritems(): # loop through tree vertices
        for qk,qv in Q.iteritems(): # loop through non-tree vertices
            dist = fv['G'].distanceTo(qv['G']) # calculate distance
            if dist < min_dist: # if distance is less than current minimum, remember
                fk_fv_qk_qv = [fk,fv['G'],qk,qv['G']]
                min_dist = dist
    F[fk_fv_qk_qv[2]] = {'G': fk_fv_qk_qv[3]} # add to tree vertices
    del Q[fk_fv_qk_qv[2]] # delete from non-tree vertices
    return arcpy.Polyline(arcpy.Array([fk_fv_qk_qv[1].centroid,fk_fv_qk_qv[3].centroid]),sr) # return new line

fc = 'points' # points feature class
sr = arcpy.Describe(fc).spatialReference # spatial ref
Q = {str(i[0]):{'G':i[1]} for i in arcpy.da.SearchCursor(fc,['OID@','SHAPE@'],spatial_reference=sr)} # non-tree vertices
F = {} # empty tree vertices
lines = [] # placeholder for lines
q_cur = Q.keys()[0] # get first non-tree vertex
F[q_cur] = Q[q_cur] # transfer to tree vertices
del Q[q_cur] # delete from non-tree vertices
while len(Q)>0: # do until all non-tree vertices assigned to tree
    output = addToTree(F,Q,sr) # add a vertex via function
    lines.append(output) # remember line
arcpy.CopyFeatures_management(lines,r'in_memory\lines') # write lines


texts: Custom renderer crashes QGIS because "renderFeature" is called in a different thread

texts:
I am trying to implement a custom renderer with Python in QGIS. I have put together a small toy example based on the code from the Creating Custom Renderers section of the PyQGIS Cookbook.

To get this code to work, I had to override the renderFeature function from the QgsFeatureRenderer class. I believe this should not be necessary, because I am essentially just re-implementing the original implementation. However, if I remove the implementation of renderFeature from my code and try to use the renderer in a version of QGIS that was built using CMAKE_BUILD_TYPE=Debug, the following assert failure is thrown:

Fatal: ASSERT failure in QgsFeatureRenderer::renderFeature: "renderFeature called in a different thread - use a cloned renderer instead", file ../src/core/symbology/qgsrenderer.cpp, line 122

Note that in a QGIS version without the Debug build type, the code works even with the renderFeature function removed. I assume the assertion is there for a reason, though, so ignoring this seems like the wrong thing to do.

What am I missing? Which renderer do I have to clone to avoid this error?

Steps to reproduce:

Start a version of QGIS built with the Debug build type
Paste the code below into the Python console and hit Enter
Load a points layer
In the layer styling panel, select the Test renderer
QGIS should crash

from qgis.core import QgsWkbTypes, QgsSymbol, QgsFeatureRenderer
from qgis.gui import QgsRendererWidget
from PyQt5.QtGui import QColor

class TestRenderer(QgsFeatureRenderer):
    def __init__(self):
        QgsFeatureRenderer.__init__(self, "TestRenderer")
        self.sym = QgsSymbol.defaultSymbol(QgsWkbTypes.geometryType(QgsWkbTypes.Point))
        self.sym.setColor(QColor('blue'))

    def symbolForFeature(self, feature):
        return self.sym

    def startRender(self, context, vlayer):
        self.sym.startRender(context)

    ####
    # The following function is basically identical to the default implementation of renderFeature
    # However, if I don't implement it here, QGIS crashes and I get the following
    # error message:
    # "renderFeature called in a different thread - use a cloned renderer instead"
    ####
    #def renderFeature(self, feature, context, layer, selected, drawVertexMarker):
    #   s = self.symbolForFeature(feature)
    #   s.renderFeature(feature, context, layer, selected, drawVertexMarker)
    #   return True

    def stopRender(self, context):
       self.sym.stopRender(context)

    def usedAttributes(self, context):
        return []

    def clone(self):
        return TestRenderer()

class TestRendererWidget(QgsRendererWidget):
    def __init__(self, layer, style, renderer):
        QgsRendererWidget.__init__(self, layer, style)
        if renderer is None or renderer.type() != "TestRenderer":
            self.r = TestRenderer()
        else:
            self.r = renderer

    def renderer(self):
        return self.r

from qgis.core import QgsRendererAbstractMetadata, QgsApplication

class TestRendererMetadata(QgsRendererAbstractMetadata):
    def __init__(self):
        QgsRendererAbstractMetadata.__init__(self, "TestRenderer", "Test renderer")

    def createRenderer(self, element):
        return TestRenderer()

    def createRendererWidget(self, layer, style, renderer):
        return TestRendererWidget(layer, style, renderer)

QgsApplication.rendererRegistry().addRenderer(TestRendererMetadata())


texts:
I would like to use Ipython in the console to perform tasks more easily.

How can I do that?

I want to use Ipython because I prefer to program with it. I would like to have TAB-completion because the attribute of QGIS module made me crazy.


texts: Finding point from Layer B that matches point from Layer A using search radius

texts:
I have two point shapefiles of trees from two different surveys. I want to verify if data from Survey A matches data from Survey B.

I want to create a search radius of 10m around each tree from Survey A, to find if there is a tree in Survey B that is of the same species (height or other parameters) within that search radius.

So far, the method I've come up with is to create a Survey A buffer before spatially joining it with Survey B points.

Is there a better way to do this?

I am using ArcMap 10.6.1.


texts:
First time using OpenLayers.  I am getting an error "goog is not defined".  Here are the steps I performed:

Downloaded and unzipped v3.0.0-gamma.4
It created the following directory structure:

apidoc
build
closure-library
css
doc
examples
ol
resources

I proceed to create a file called ol3.html and place it in the examples folder shown above.  Here is the code of that file:

<!doctype html>

    <link rel="stylesheet" href="../css/ol.css" type="text/css">

    <title>OpenLayers 3 example</title>

My Map

When I load up the page, I get the error "goog is not defined".  I can see that it finds the local ol.js file but it give the error "goog is not defined".
If I switch this line:

To this:

Then everything is fine and I can see the map. What am I missing?


texts: GeoWebCache Layer name is appended with SRS

texts:
I am trying to display ESRI Tiled Cache on the browser using GeoWebCache and OpenLayers. I am following these instructions to serve arcgis exploded cached tiles as a service.I have successfully configured all the xml, and the layer is showing on GeoWebCache Portal, but it is also adding layer name in SRS of geowebcache layer. Here is screenshot of Geowebcache layers

This is my XML configuration added in geowebcache.

<arcgisLayer>
      <tilingScheme>D:\\path_to_cache\\conf.xml</tilingScheme>
        <name>fsd</name>
</arcgisLayer>

Due to this, i am not able to see the layer in browser if i view it in by clicking on .png open layer URL.


texts: What is the easiest one click install method for PostGIS and it's entire dependencies?

texts:
For someone new to database servers and their dependencies, what is the easiest assured way of installing PostGIS on a machine with no other software beyond what comes with its operating system? I'd like to know for Ubuntu, Windows and if possible a common Linux method. I'm understanding that you need to install PostGreSQL first, but before that perhaps you need to install Apache Server, and I'm not sure if that has dependencies.

Perhaps there a single binary that loads EVERYTHING for Windows, and maybe there a single program from the Ubuntu GIS repositories that installs everything. If not, what's the install order?


texts:
I am trying to get a map with the type of soil depending of their texture. I have tried three options:

In SAGA.
By using the tool "soil texture classification" where by using three raster with the % of sand, silt and clay it was supposed that I was going to get a map with results between 1 - 12 (http://www.saga-gis.org/saga_tool_doc/2.2.1/grid_analysis_14.html), but what I get it is a map with data between 254 and 255. Nothing to do with the expected result.

In R.
There is a package named "soiltexture", but when I try to install it I am giving this message: "ERROR: dependency ‘sp’ is not available for package ‘soiltexture’". I do not know how to fix this issue.

In GRASS.
It seems that there is a tool named "r.soils.texture", but I do not have this tool in my grass (QGIS 2.18.15).

So, I do not see more options to be able to get a map of soils by using any of these software.

Has anyone had any of these problems before?


texts:
I have a Landsat footprint polygon shapefile that I simply want to create a 6 x 6 fishnet on. I thought this would be simple, but every time I run the tool I get a fishnet output that extends longer than the original footprint and is shorter in height. I think the problem may stem from the fact that the footprint is slanted and the y-axis may be a different angle than the x-axis.

I used Feature Vertices to Points and then Add XY Coordinates to get the exact coordinates of the footprint polygon's corners. Then within the Fishnet tool, I used the footprint's extent as the 'Template Extent', the bottom left corner as the 'Origin Coordinate', the top left corner as the 'Y-Axis Coordinate', and the top right corner as the 'Opposite corner of Fishnet.'

Coordinate system is WGS_1984_UTM_Zone_18S.

Any ideas for how I can fix this and generate a properly placed fishnet?

Here's the strange looking fishnet output:

And here are the parameters I passed into the fishnet tool, along with the proper corner coordinates in the top left of the image (Point_X and Point_Y):


texts: Why doesn't PostGIS use gist index when doing a ST_DUMP(ST_UNION)?

texts:
My query:

DROP TABLE IF EXISTS tmp;
CREATE TEMP TABLE tmp AS SELECT *, ST_BUFFER(the_geom::GEOGRAPHY, 3000)::GEOMETRY AS buffer FROM af_modis_master LIMIT 20000;
CREATE INDEX idx_tmp_the_geom ON tmp USING gist(buffer);
EXPLAIN SELECT (DUMP(ST_UNION(buffer))).path[1], (DUMP(ST_UNION(buffer))).geom FROM tmp;

Output from EXPLAIN:

Aggregate  (cost=1705.52..1705.54 rows=1 width=32)
  ->  Seq Scan on tmp  (cost=0.00..1625.01 rows=16101 width=32)

Seq Scan means it is not using the index, right? Why not?

UPDATE: Even adding a where clause that filters based on the buffer causes a Seq Scan:

ANALYZE tmp;
EXPLAIN SELECT (DUMP(ST_UNION(buffer))).path[1], (DUMP(ST_UNION(buffer))).geom FROM tmp WHERE ST_XMIN(buffer) = 0.0;

(This question was also posted here. Apologies for cross-posting but the community there is a bit more active (maybe just on PostGIS query topics?), so perhaps wil provide an answer quicker.)


texts: My OSM country roads network is not routable

texts:
I need a routable shapefile in order to execute pgrouting algorithms. I download a national zip from http://download.geofabrik.de/europe/greece-latest.shp.zip and import the roads shapefile in my PostgreSQL/PostGIS database through QGIS. Then I execute the following SQL commands (following the workshop instructions):

ALTER TABLE roads ADD COLUMN length integer;
UPDATE roads SET length = st_length(st_transform(geom,2100));
ALTER TABLE roads ADD COLUMN reverse_cost double precision;
UPDATE roads SET reverse_cost = length;
ALTER TABLE roads ADD COLUMN "source" integer;
ALTER TABLE roads ADD COLUMN "target" integer;
SELECT pgr_createTopology('roads', 0.0005, 'geom', 'id');
SELECT seq, id1 AS node, id2 AS edge, cost FROM pgr_dijkstra('
            SELECT id AS id,
                     source::integer,
                     target::integer,
                     length::double precision AS cost
                    FROM roads',
            number1, number2, false, false);

where number1 and number 2 are values I define manually. All commands are executed with success but the last one. The query does not work in most cases, although there is an apparent road link, because, as I zoom in the layer, I can see breaks in the network.
I searched similar problems in the forum and I found suggestions that I should fix my shapefile in order to make it routable. I tried the GRASS commands split and clean but I met problems. With split I had no result since it said it found too many errors. With clean I had a cleaned shapefile but a non-routable one again.

Is there another way to make my dataset routable?


texts:
I'm looking for a database that contains all (most, many...as many as it can) world airports.  I know of one such set available - http://www.world-airport-database.com/database.php - but the terms state that it's for internal use only...which kind of defeats the purpose.

The data would need to have the airport code, name, location (country, city, state/locale, zip/postal, lat/lon).

Do you know of such a set?  Certainly willing to purchase such a set, free is always great too.


texts: How to trace SQL queries sent by ArcGIS Server (ArcSDE) to Oracle database?

texts:
I'd like to generate a log file containing all SQL queries sent by ArcGIS Server (ArcSDE) to Oracle database.
Is there a way to do so?
I'm using Oracle 11g and ArcGIS Server 10.0 on Windows. ArcSDE is used in direct connection.


texts: Calculate area intersected areas, problem

texts:
(I'm working in QuantamGIS)

I intersected two layers creating a new layer. Within the attibute table of this new intersect layer i would like to calculate the new areas.

However, it just doesn't calculate the correct area size. The layers are all in the same CRS and the canvas units are in meters. However I am unable to recalculate the correct area.

To check if the problem was created by me I created a polygonal area, calculated the size and came out to about the right size, however the same polygon area within the intersect layer just gives really inconsistent small numbers.

Is there any good tutorial or procedure that describes the intersect action, thus creating new polygons, with the right CRS and then recalculating the areas?


texts: Why can't I download the georeferencing plugin for QGIS 1.7.4?

texts: I used the georeferencing plugin with QGIS 1.6 and now that I've upgraded to 1.7.4 I'd like to do it again. But when I search for the georeferencer plugin, it's not available! Am I doing something wrong?... I also seems that there's no georeferencing option included with version 1.7.4.


texts: Database modeling for big data (big set of points X big sets of possible values for each points)?

texts:
I'm building a PostgreSQL database with the intent to host modelisation results for different projects.

Each project has a set of receptors points (~5000-10000 receptors) and have many results sets (could be as much as 100 000 results sets (here called scenarios results)).

Next figure is my actual (simplified) design:

The problem with this design is that the table receptor_results is going to be extremely heavy. For a project with 5000 receptors and 50 000 scenarios, it will generate 250 000 000 rows!

Any suggestion on how I can avoid such thing ?


texts: Joining point data with national census-tract layer in ArcMap

texts:
I have downloaded national census-tract level data, specifically Walkability Index, to determine within what category of this variable my point data (i.e. my sample) are classified. I have opened the layer in ArcMap without a problem, but cannot figure out how to join this layer with my "Final Sample" data. When I right click the Walkability layer, there is no attribute table or option to join data.
How do I go about this analysis if these layers cannot interact with each other?
I am very close to zooming in on each of almost 500 points in my sample and coding by hand.


texts: R - cartography package - position of legend?

texts:
I'm using the cartography package in R to plot a choropleth map. However, I have a problem placing the legend, which can be placed by means of the legend.pos argument, but which always overlaps with my map. The documentation says:

  legend.pos    : position of the legend, one of "topleft", "top",
  "topright", "right", "bottomright", "bottom", "bottomleft", "left" or
  a vector of two coordinates in map units (c(x, y)). If legend.pos is
  "n" then the legend is not plotted.

(Source: https://cran.r-project.org/web/packages/cartography/cartography.pdf)

Using "right", "bottomright" etc, the legend still overlaps my map; thus, I'd like to use the coordinates version. However, I tried already lots of different combinations, but the legend does not appear in the map anymore.

My question is: How to understand "coordinates in map units"? How can I find the right coordinates to that the legend becomes visible?

Here is my code; as you can see, I use the standard "nuts0" map.

choroLayer(spdf = nuts0.spdf, df = nuts0.df, var = "cdp",
           breaks=seq(from=0, to=240, by=40),
           col = farb,
           colNA="grey70",
           border = "grey40", lwd = 0.5,
           legend.pos = "bottomright",
           legend.title.txt = "No. companies",
           add = F)


texts:
Is it normal that geemap doesn't display attributions in my notebook? Is there a way to add them?
I'm running :
import geemap

m = geemap.Map(center=[40,-100], zoom=4)
m.add_basemap('CartoDB.Positron')
m

And I get this map


texts: Measuring distance in ArcGlobe using Python AddIn

texts:
I want to make a Python script that will calculate the distance of the line I drew and output it in a special form, but I can't figure out what I'm doing wrong.
import arcpy
import pythonaddins

class ToolClass2(object):
    """Implementation for probn_addin.tool (Tool)"""
    def __init__(self):
        self.enabled = True
        self.shape = "Line" # Can set to "Line", "Circle" or "Rectangle" for interactive shape drawing and to activate the onLine/Polygon/Circle event sinks.
    def onMouseDown(self, x, y, button, shift):
        pass
    def onMouseDownMap(self, x, y, button, shift):
        print(x,y)
        return(x,y)
    def onMouseUp(self, x, y, button, shift):
        pass
    def onMouseUpMap(self, x, y, button, shift):
        print(x,y)
        return(x,y)
    def onMouseMove(self, x, y, button, shift):
        pass
    def onMouseMoveMap(self, x, y, button, shift):
        pass
    def onDblClick(self):
        pass
    def onKeyDown(self, keycode, shift):
        pass
    def onKeyUp(self, keycode, shift):
        pass
    def deactivate(self):
        pass
    def onCircle(self, circle_geometry):
        pass
    def onLine(self, line_geometry):
        onMouseDownMap(x1,y1)
        onMouseUpMap(x2,y2)
        dist = math.sqrt((x2 - x1)**2 + (y2 - y1)**2)
        message = "Your lenght:" + str(dist)
        pythonaddins.MessageBox(message, "Ruler")
    def onRectangle(self, rectangle_geometry):
        pass


texts: ArcMap 10; Syntax for Minimum operator in raster calculator

texts:
I want to linear stretch my raster grid formate having column/rows 3128, 3503 with 30 meter resolution by using the following syntax in raster calculator of arcmap 10.1.

("raster" - FocalStatistics("raster", NbrRectangle, 30, 30, "Minimum") /
(FocalStatistics("raster", NbrRectangle, 30, 30, "Maximum") -
FocalStatistics("raster", NbrRectangle, 30, 30, "Minimum"))

But this gives me the error while executing.
I have also tried the simple equation but did not work.

("raster" - min("raster")) * 1.0 / (max("raster") - min("raster")) + 0.0

Error displayed:

  000539 : Error message from Python.
  ArcGIS 10.1 Locate topic
  Description
  The calculation used by the Calculate Field or Calculate Value tool is
  invalid. This error message provided will list the specific Python
  error.
  Solution
  This error code covers a number of Python errors:
  Example error 1: exceptions.TypeError: cannot concatenate 'str' and 'int' objects.The
  above is a Python-specific error. The calculation is attempting to add
  or concatenate a string and a number.
  Example error 2: Invalid field shape@distance
  The above is an error using the geometry object. The distance method
  is not a valid method of the geometry object.For specific Python
  issues, consult the external Python help for more information, or
  consult the Calculate Field or Calculate Value help for more
  information on these tools.

What is the right syntax.


texts:
I'm collecting POIs data from several sources and plotting them on map, i'm having difficulties merging the POIs because some of them are redundant between the data sources i'm using. How can i assign a unique identifier for each POI in a ways that the implementation logic can figure out the redundancy, e.g. I have a starbucks in two of the sources each for the same branch, and there is slight difference in the long/lat value due to accuracy, how can i automatically identify th redundancy.

Street names are not enough as they are not unique and might get updated from different sources, i was thinking of using the lat/long in addition to the first three letters of the POI, something like STA-31,444-30,122 but again the long/lat are not always the same between the different sources because of the accuracy, they might be off for a few meters.

Thanks
Yehia


texts:
I'm trying to generate a GeoJSON  with this code :
 @RequestMapping(value = "/testpoly", method = RequestMethod.GET)
  @ResponseBody // tells a controller that the object returned is automatically serialized into
                // JSON and passed back into the HttpResponse object.
  public JSONObject CoucheTest() {

    List<String> ports = repositorycouche.selectAll();

    JSONObject featureCollection = new JSONObject();
    featureCollection.put("type", "FeatureCollection");
    JSONObject properties = new JSONObject();
    properties.put("name", "ESPG:4326");
    JSONObject crs = new JSONObject();
    crs.put("type", "name");
    crs.put("properties", properties);
    featureCollection.put("crs", crs);
    JSONArray features = new JSONArray();

    for (int i = 0; i < ports.size(); i++) {

      JsonObject port = (JsonObject) new JsonParser().parse(ports.get(i));
      JSONObject feature = new JSONObject();
      feature.put("type", "Feature");
      feature.put("geometry", port);
      features.put(feature);

    }

    featureCollection.put("features", features);

    return featureCollection;
  }

}

It's giving me as a result :
{"map":{"features":{"myArrayList":[{"map":{"geometry":{"type":"MultiPolygon","coordinates":[[[[-5.79241275787353,35.792092328535],[-5.80159664154053,35.7995762437594],[-5.82146644592285,35.7915005470267],[-5.80185413360596,35.7745458760594],[-5.78661918640137,35.7864180067475],[-5.79241275787353,35.792092328535]]]]},"type":"Feature"}},{"map":{"geometry":{"type":"MultiPolygon","coordinates":[[[[-5.50157546997069,35.9122717133218],[-5.60337066650389,35.8466217303553],[-5.57865142822264,35.8215713916119],[-5.5462074279785,35.8332625339134],[-5.46432495117186,35.8872421214509],[-5.50157546997069,35.9122717133218]]]]},"type":"Feature"}}]},"crs":{"map":{"type":"name","properties":{"map":{"name":"ESPG:4326"}}}},"type":"FeatureCollection"}}

How can I remove "map" and "myArrayList" from my result?
Ps : I'm using org.json package.


texts: QGIS script doesn't commit attribute changes

texts:
I have a script to generate random points within a polygon, and to give the points the attributes of their host polygon. It generates points, but for some reason it won't commit the changes to the attributes. The weird thing is that the output while running the script clearly shows the feature it is adding has the attributes, but when it's done I'm left with a bunch of featureless points with only an ID.

I have a feeling it's something fundamental I'm missing in committing the attributes.

import random

# Prepare new temporary editable memory layer
pointLayer = iface.addVectorLayer("Point?crs="+source.crs().toWkt(), "random_points", "memory")

for f in source.getFeatures():
    # Create bounding box from polygon
    bounds = f.geometry().boundingBox()
    attributes = f.attributes()
    xmin = bounds.xMinimum()
    xmax = bounds.xMaximum()
    ymin = bounds.yMinimum()
    ymax = bounds.yMaximum()

    for i in range(['Number']):

        inBounds = False
        while not inBounds:
            # generate random point 'p' within feature bounds
            xRandom = xmin + (random.random() * (xmax-xmin))
            yRandom = ymin + (random.random() * (ymax-ymin))
            p = QgsPoint(xRandom,yRandom)
            pGeom = QgsGeometry.fromPoint(p)

            # check if new point is contained in poly
            # if within bounds, add attributes and exit loop
            if f.geometry().contains(pGeom):
                pointFeature = QgsFeature()
                pointFeature.setGeometry(pGeom)
                print('source atts: ')
                print(attributes)
                pointFeature.setAttributes(attributes)
                print('point atts: ')
                print(pointFeature.attributes())
                pointLayer.dataProvider().addFeatures([pointFeature])
                inBounds = True


texts: How to Update attribute values from queried features on SimpleFeatureCollection?

texts:
We want to update some features of a Layer. We have used the following Code to update features but not successful.

Firstly Get SimpleFeatureSource by a Feature Layer.

SimpleFeatureSource sfSource = (SimpleFeatureSource) Lyr.getFeatureSource();

After quried on it we got some data in the form of SimpleFeatureCollection.

org.opengis.filter.Filter flt = ff.within(ff.property("the_geom"),ff.literal(polygon));
SimpleFeatureCollection sfc = sfSource.getFeatures(flt);
SimpleFeatureIterator iterator = sfc.getFeatures().features();
while (iterator.hasNext()) {
   SimpleFeature feature = iterator.next();
   feature.setAttribute(0, object);
}

If we do this feature's attribute will update and reflect in sfc (SimpleFeatureCollection) but not reflect in sfSource (SimpleFeatureSource).

Updated:-

Transaction transaction = new DefaultTransaction("Transaction");
    Maplayer MpLayer = this._MapListLayers.get(cmbLayers.getSelectedItem().toString());
    FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2();
    Filter flt = ff.id(ff.featureId("long cl0805143437.44177"), ff.featureId("long cl0805143437.44178"),
            ff.featureId("long cl0805143437.44179"), ff.featureId("long cl0805143437.44180"),
            ff.featureId("long cl0805143437.44181"), ff.featureId("long cl0805143437.44182"),
            ff.featureId("long cl0805143437.44183"), ff.featureId("long cl0805143437.44184"),
            ff.featureId("long cl0805143437.44185"), ff.featureId("long cl0805143437.44186"),
            ff.featureId("long cl0805143437.44187"), ff.featureId("long cl0805143437.44188"),
            ff.featureId("long cl0805143437.44189"), ff.featureId("long cl0805143437.44190"));
    try {
        MpLayer._SimpleStore.setTransaction(transaction);
        MpLayer._SimpleStore.modifyFeatures("Equipment", new String("Equipment"), flt);
        transaction.commit();
    } catch (Exception e) {
        e.printStackTrace();
        transaction.rollback();
    }

As per your comments we have tried upper code but results me in an error shown below:-

  java.lang.RuntimeException: OGR Error: General Error
      at org.gdal.ogr.ogrJNI.Layer_SetFeature(Native Method)
      at org.gdal.ogr.Layer.SetFeature(Layer.java:145)
      at OgrSource.JniOGR.LayerSetFeature(JniOGR.java:284)
      at OgrSource.OGRDirectFeatureWriter.write(OGRDirectFeatureWriter.java:114)
      at org.geotools.data.InProcessLockingManager$1.write(InProcessLockingManager.java:337)
      at org.geotools.data.store.DiffTransactionState.commit(DiffTransactionState.java:161)


texts: Raster and digitized shapefile misplaced in same projection

texts: I have georeferenced .tiff rasters with cadastral borders and corresponding digitized borders in shapefile based on that raster. All that was done in some Mapsoft software. When I insert both in QGIS in same projection WGS84/UTM, they are slightly misplaced for some 0.5 m in NW-SE direction, same is in ArcMap. Any hint why that hapenning?


texts: Geoplot/Cartopy Orthographic projection won't work near equator

texts:
I am having trouble getting Geoplot's Orthographic projection to work near the equator.  It works fine at locations ±90 longitude where latitude is right at 0, but if I run the code below, it fails the try blocks if -90 < longitude < 90 when latitude is at or near 0.  In addition, at the prime meridian, the try blocks will fail if latitude is as high as 6.5 or as low as -0.2.

Why?  What's wrong with the projection?

The problem with Geoplot is the same when I use Cartopy directly, as Geoplot is built on top of Cartopy.

#Imported packages for the geographic part of the program
import geopandas as gpd
import matplotlib.pyplot as plt
import urllib.request as ur
import geoplot as gpl
import geoplot.crs as gcrs
from shapely.geometry import Point

#Imported packages for the graphical interface part of the program
import tkinter as tk
from tkinter import *

#Function for the geographic part of the program
def show_map(lon, lat):
    #Parameters (longitude and latitude centres of the projection)
    lonctr = lon
    latctr = lat
    centrepoint = Point(lonctr, latctr)

#<snip>----------------------------------

    #Creating a world map and ensuring the map is limited by the orthographic
    #projection
    themap = plt.figure("World Map")

    #Preparing projected data for the Natural Earth file to be plotted
    #on the axes using geoplot
    try:
        check = 1
        themap = plt.subplot(111,projection=gcrs.Orthographic(central_longitude=lonctr,
            central_latitude=latctr))
        land = gpl.polyplot(earth, projection=gcrs.Orthographic(central_longitude=lonctr,
            central_latitude=latctr), linewidth=1, edgecolor="white", facecolor="blue",
            zorder=1, ax=themap)
        dot = gpl.pointplot(newdata, projection=gcrs.Orthographic(central_longitude=lonctr,
            central_latitude=latctr), s=20, color="red", marker="o", zorder=4, ax=themap)
    except:
        themap = plt.subplot(111,projection=gcrs.PlateCarree(central_longitude=lonctr))
        land = gpl.polyplot(earth, projection=gcrs.PlateCarree(central_longitude=lonctr),
            linewidth=1, edgecolor="white", facecolor="blue", zorder=1, ax=themap)
        dot = gpl.pointplot(newdata, projection=gcrs.PlateCarree(central_longitude=lonctr),
            s=20, color="red", marker="o", zorder=4, ax=themap)

    #Using projected data from overlay file to be plotted on the axes
    #using geoplot
    try:
        if(check == 1):
            polyover = gpl.polyplot(overlay, projection=gcrs.Orthographic(central_longitude=lonctr,
                central_latitude=latctr), linewidth=1, edgecolor="white", facecolor="pink",
                zorder=2, ax=themap)
        else:
            polyover = gpl.polyplot(overlay, projection=gcrs.PlateCarre(central_longitude=lonctr), linewidth=1,
                edgecolor="white", facecolor="pink", zorder=2, ax=themap)
    except:
        pass
    try:
        if(check == 1):
            pointover = gpl.pointplot(overlay, projection=gcrs.Orthographic(central_longitude=lonctr,
                central_latitude=latctr), s=20, color="green", marker="o", zorder=3, ax=themap)
        else:
            pointover = gpl.pointplot(overlay, projection=gcrs.PlateCarree(central_longitude=lonctr),
                s=20, color="green", marker="o", zorder=3, ax=themap)
    except:
        pass

    check = 0

#<snip>----------------------------------


texts: How to make offsetting points perpendicular to a line (without ET GeoWizards)?

texts:
So I'm trying to create offsetting points perpendicular from a line and have run into a roadblock in Arc 10.1. I have the ability to use the ET GeoWizards tool 'Points Along Polylines' but the tool only functions correctly in the ET GeoWizards window (not in toolbox, ModelBuilder or Python which is where I want to run it from). Essentially, I would like to have the same results has the ET GeoWizards tool 'Points Along Polylines'.

I had tried experimenting with creating a set buffer around my Polyline features and using the  several different tools to try to snap a line to the extent of the buffer...but some of my line features have drastic angles which creates a buffer overlap and doesn't produce the results I want.

The final product is within a Python script, so outside of the box thinking is appreciated!


texts: Convert Layer to Shapefile without losing Feature definitions

texts:
I have created a worldmap with custom defined features, representing State/Zip-Code sales areas. I want to convert them into a Shapefile, so I can convert the Shapefile into a polygon definition file for usage in Tableau.

When I do the following conversion, I lose the Feature definitions:
- I select the layer I want to export
- I do a DXF export
- I run Dxf2Shp Converter (plugin) to convert the DXF to a SHP
- When I then convert the SHP file to a Polygon definition file, I lose the feature groupings/definition.

A different Polygon shape is created for each little Island. For example, Hawaii would be separated into many different individual polygon shapes, while in my QGIS features, they are all grouped together as one. Basically, I'm losing the Feature grouping and get different polygons for each geographical island.

Any idea how to change/correct this? Or is there a better way of converting a QGIS file into a Shapefile?


texts: Adding North Arrow in ArcGIS Online web application/map?

texts:
I would like to add a north arrow for online applications/maps in the ArcGIS Online platform.

I can't seem to find if this is possible and/or how?


texts: Choosing coordinate system to store geography data for celestial coordinates

texts:
I am doing an astronomy project. I want to have the information about our images stored in a spatially enabled database. This, I would think, should be a very easy special case for GIS functions because the sky can be treated as perfectly spherical, and does not require an elliptical treatment like the surface of the earth. Unfortunately, I have not found a way to do this yet and I have been dodging mines with spatial functions that use an elliptical earth. (Pretty much any function that returns meters instead of degrees may be using an elliptical calculation. Luckily, many of the PostGIS functions I've needed appear to have incomplete implementations where the documentation explicitly states that the returned results are for a sphere and not for the ellipsoid. But that may change with future versions, which is a cause for concern.)
Background: I am presently using PostgreSQL with PostGIS and WGS 84 coordinates (SRID=4326). This works fairly well. I am creating a closed POLYGON from the right ascension and declination of the four corners of the image. I have a lot of images (10k or more), covering a large area of the sky. Each images is about 1 degree square. From the set of these images, I am making mosaics from small subsets of 15 to 30 images. Each mosaic is about 1.5 degree square.
Presently, I am storing the geography of the mosaics as a MULTIPOLYGON that consists of all the POLYGONS corresponding to each image that went into the mosaic. [A better solution would be to create a single POLYGON that describes the perimeter of the union of all the individual polygons. I don't know if this can be done in spherical coordinates (i.e., that the geography type). This would also be an interesting answer for me too.] The date line and celestial poles may be included in an image in the dataset so I have been avoiding projecting to planar coordinates to the extent possible.
What coordinate system ought I use for celestial coordinates with PostGIS functions?
I have looked at  http://spatialreference.org/ but have not found anything so far. Google has turned up little. I am stumped. Basically, I want to ensure that if a function returns meters as a distance, it is meters along a great circle on a sphere.

I am using PostGIS 1.5.2. I have not yet tried PostGIS 2.0. I am curious if the ST_CoveredBy function works with a POLYGON and a MULTIPOLYGON of type geography. If anybody is running 2.0, could you tell me if you get the same error as this:
mydb=# select ST_CoveredBy(ST_GeographyFromText('MULTIPOLYGON(( (10.37795 -69.57926,8.9498 -69.54875,9.0178 -69.21643,10.4242 -69.24648,10.37795 -69.57926),(10.42436 -69.24618,9.01774 -69.2162,     9.08363 -68.88389,10.46914 -68.91344,10.42436 -69.24618)))'),ST_GeographyFromText('POLYGON((10.46915 -68.91315,9.08371 -68.88364,9.14755 -68.5513,10.5125 -68.58038,10.46915 -68.91315))'));
ERROR:  geography_covers: only POLYGON and POINT types are currently supported
CONTEXT:  SQL function "st_coveredby" statement 1

I have tried PostGIS 2.0. This function still only works on points and polygons, not more general shapes.


texts: Using ArcMap, how may I convert a shapefile to use decimal degrees?

texts:
I have a shapefile that currently uses meters as the measure for its coordinate system, but the map framework I wish to use expects decimal degrees. I'd like to reproject my shapefile to use decimal degrees, but I can't seem to figure out how to do so.

In ArcMap, I bring up the Project tool (ArcToolbox -> Data Management Tools -> Projections and Transformation -> Feature -> Project), and I select my desired projection (Projected Coordinate Systems -> Continental -> North America -> USA Contiguous Albers Equal Area Conic USGS.prj), but I'm not sure how to select my desired units at that point. If I click "Modify..." I can edit the Projected Coordinate System Properties, and I would expect "decimal degrees" to be an option in the "Linear Unit" dropdown, but it's not.

How may I go about choosing to reproject into decimal degrees?

EDIT:

Short answer: I wasn't understanding the difference between Geographic Coordinate Systems and Projected Coordinate Systems. A Projected system describes a flat surface, and so it doesn't make sense to describe it in terms of degrees; whereas a Geographic system uses a 3-dimensional spherical surface, so degrees are appropriate.

A system like the Google Maps, which uses the Web Mercator projection along with latitudes and longitudes, presumably must reproject lat/lng coordinates on the fly.

Thanks to @Brad Nesom for helping to elucidate me. I'm accepting his answer for being patient and helpful :-)


texts: Generating new point file based on specified distances and azimuths/angles from previous point in ArcMap?

texts:
I have a shapefile of multiple points (upwards of 100). I would like to generate a new shapefile with points that are a specified distance and direction away from the previous points. For example, each point has a unique distance and angle/bearing associated with it. I have read into the create features and using the distance/direction tool. However, this is very tedious as I have many points, and will continue to. I have also found the "Bearing Distance to Line" tool and I believe this does exactly what I want it to. however, it generates a line, and not a point.

Is there a way to do it that makes a line, or a way to make new shapefile with only points at the end of each line?


texts:
Trying to add own graphics (.svg file containing lines as an extra layer) to the map in QGIS.

How to do it?


texts: TileMill: Can't use both background-color and background-image?

texts:
Very new to TileMill and am working on my first project. I'm using a pattern image to enhance some areas on the map, water being an example.  I've added the pattern to the water areas on my layers, and am now trying to match the background to those areas.  Having trouble - can someone please tell me what I'm doing wrong?

This renders only the pattern in the map preview, but not the background color as well:

Map {
  background-color:@Water;
  background-image:url(images/water.png);
}


texts:
I am having a question about the raster calculation as I am trying to create a tool which I want to link with my python script to sum weight rasters but if failed to do so. here is my code:
     #*the original code is very long so I will try to minimize to my questions  *

import arcpy
arcpy.env.workspace = "c:my workspace"
from arcpy.sa import *
...
input_weights = arcpy.GetParameterAsText(0)  ***#this is my weights: in model builder, the DATATYPE is the DOUBLE #***
input_weights2 = arcpy.GetParameterAsText(1)***#this is my weights:in model builder, the DATATYPE is the DOUBLE #***
RASTER_1 =arcpy.GetParameterAsText(3) #my feature class
RASTER_2 =arcpy.GetParameterAsText(4) # my feature class
...

myraster  = EcuDistance(RASTER_1)
myraster2 = EcuDistance(RASTER_2)
weight =myraster*input_weights+myraster2*input_weights2
weight.save("xxx")
......

Every time the tool would stop at (weight =myraster*input_weights+myraster2*input_weights2) and it says error:can't multiply sequence by non-int of type 'float'?

so if I manually entered those weights like:

input_weights = .5
input_weights2 = .4
...
...
weight =myraster*.5+myraster2*.4
weight.save("xxx")

the tool worked perfectly.Is there something wrong with my code. Should I do a transformation for the weights from input of the tool


texts: How does manual Label Placement in QGIS 1.9 work?

texts:
There are a bunch of questions here related to labeling in QGIS, but it's not clear which version they are referring to, and I know this is something that has changed recently. Hopefully it makes sense to start a new question due to the new labeling functionality.

I'm trying to label a road network so that every street is labelled. I have the 'auto' labeling doing a decent job, but it is not placing a label on every street. How can I (with the minimum of effort) keep the existing auto-placed labels (possibly with a few tweaks/manual moves,) but then also go in and manual 'force' missing labels to appear, and position them myself and/or change font size to make them fit?

The new label tools in the toolbar are not very easy to figure out, and the 'pinning' concept is lost on me. Can anyone explain the steps?


texts: ArcGIS Server Image Service from Mosaic Layer

texts:
I have successfully created an image service from a mosaic dataset using NAIP imagery. I would like to create another image service that displays the same area in NIR (4,3,2). I have created a mosaic layer with the display using 4,3,2 and I was able to use the layer to create a new image service.(pointing to the same mosaic dataset) But the new image service does not display in NIR. I have also tried generating tile cache from the NIR mosaic layer and copying to the server, and it still does not display.

If I add the mosaic layer into arcmap, and open the attribute table for the mosaic dataset and select to view individual rasters, they do display in NIR.

The point of me doing this is so I can display my imagery in different ways, while referencing the same image files, to save space.

Has anyone successfully accomplished this?


texts: Get feature values from ImageCollection to list

texts:
I am calculating the mean precipitation value over a given area from CHIRPS precipitation data.
My goal is to get a list containing a date and precipitation value for every Image of the ImageCollection.
Using .getInfo() on the reduced ImageCollection I can get a dict representing the ImageCollection and from there select the values on the client side.
import ee
from datetime import datetime

ee.Initialize()

dataset = ee.ImageCollection('UCSB-CHG/CHIRPS/DAILY').filter(ee.Filter.date('2016-01-01','2016-12-31'))
area = ee.Geometry.Polygon([[[29.045341997822412, -2.1339228457039368],
                                               [29.07984593458999, -2.1339228457039368],
                                               [29.07984593458999, -2.0966977113182073],
                                               [29.045341997822412, -2.0966977113182073],
                                               [29.045341997822412, -2.1339228457039368]]])

 # Image reduction applied to each image.
def reduce_dataset_region(image):
    # Calculate mean of precipitation on defined area.
    local_precipitation_image = image.reduceRegion(
        reducer=ee.Reducer.mean(),
        geometry=area,
        scale=20
    )

    return image.set('mean', local_precipitation_image)

# Apply region reduction to ImageCollection
reduced_dataset = dataset.map(reduce_dataset_region, True)

# Request server-side ImageCollection as dict
reduced_dataset_dict = reduced_dataset.getInfo()

# Create list to hold daily precipitation values.
precipitation_list = []

# Loop through features in image collection dict. One feature represents one image.
for feature in reduced_dataset_dict['features']:
    # Get date and turn into datetime
    # Get precipitation value and add with datetime to list.
    precipitation_list.append((timestamp, feature['properties']['mean']['precipitation']))

print(precipitation_list)

This seems slow and not the right way to go, especially since I want to use data from multiple years.
How can I create a list containing the precipitation value for each Image and send only that, and not the whole ImageCollection, to the client? What is the better alternative to .getInfo here?


texts: Assigning values of points outside a buffer to points within a buffer

texts:
I have 10m buffers around line features and a point shapefile which has some points intersecting the buffer.
I would like to take each point that falls within the buffer and assign it the value of the point it is closest to outside of the buffer.

Anybody know how to accomplish this? I'm using version 10.1.


texts: I'm having a bit of an issue with the proportional circle plugin for QGIS. I want to be able to adjust the start size of the circle but for the life of me I can not figure it out but I was able to do it on a previous map that I worked on. I'm on a Mac using QGIS 2.8.4. I also need to label the circle with a numerical value which I have been about to do but I need to reduce the starting size of the circle.


texts: Rasterizing z values from vector grid points as centroid of cells?

texts:
I have read some Q&As related to this:

Creating raster or vector raster from DEM point file?,

Does a value in srtm dem file refer to the center or corner of a cell?

However, none specifically answer my question.

I have a number of 5x5 m grids consisting of X,Y,Z coordinates (DEM) each representing one tile, once rasterized and each 5x5 km in size (a lot of data!).
If I use gdal_rasterize in QGIS, the calculated raster contains cells between four grid points, resulting in no data gaps between the tiles. What I want is a raster that assigns the z value of the grid points as centroids of the cells. That way, I expect the no data line to disappear between the tiles.
As I will have to convert approx. 90 grid files to raster tiles, a batch process would be the best solution.

Which process to rasterize a grid can I use in QGIS to assign z values as centroids of the raster cells?


texts: How to copy a line feature into a polygon personal geodatabase feature class?

texts: I have a line shapefile with a proposed hiking trail.  I am trying to copy this feature into a larger geodatabase feature class.  The feature class is set to polygons.  Is there a way I can just convert a line to a polygon? or some way to get it into this polygon feature class.  It is just a single line, not a grid or parcel lines.  I'd like it to be a polygon that just looks like a line, just a few pixels wide.


texts: ArcGIS Pro Model Builder: inline substitution to create raster filenames

texts:
I am using Model Builder in ArcGIS Pro to set up an iterative work flow to analyse satellite images.

The input tiff files have long files names such as: 5023109_2010-06-21_RE1_3A_Analytic_clip.tif

I am trying to use inline substitution to give output rasters shorter filenames. For example:
"R"+"%Name%"[10:17], where Name is the variable representing the name of the input tiff.

I have read other discussions on this, but haven’t managed to solve the problem.


texts: Raster calculation with multiple rasters using ArcGIS Desktop?

texts:
I have three rasters:

Raster1: resolution 25x25, floating point. Value of "1" when it is true, value of -9999 when it is false (only this two values)
Raster2: resolution 25x25, floating point. Only integer values from 194 to 1738
Raster3: resolution 25x25, floating point. Decimal values are stored too, from -2544,81 to 941,66

All the rasters have the same CRS and are in TIFF format.

What I want to do is to create a fourth raster like this:

If Raster1 =1, then Raster2 (true statement), else Raster3 (false statement)

I tried this in ArcGis 10.3 Raster Calculator:

Con("Raster1" == 1,"Raster2","Raster3")

But I get this error:

Con(Raster(r"Raster1") == 1,Raster(r"Raster2"),Raster(r"Raster3"))
ERROR 000539: Error running expression: rcexec()
Traceback (most recent call last):
File "<expression>", line 1, in <module>
File "<string>", line 8, in rcexec
RuntimeError: ERROR 010240: Could not save raster dataset to F:\folder\step13 with output format GRID.

Failed to execute (RasterCalculator).

I tried the same thing using Python:

import arcpy
from arcpy import env
from arcpy.sa import *
env.workspace = "F:/folder"
Ras1 = Raster("Raster1")
Ras2 = Raster("Raster2")
Ras3 = Raster("Raster3")
outCon = Con((Ras1==1),Ras2, Ras3)
outCon.save = ("F:/folder/out")

But I get this error:

RuntimeError: ERROR 000875: Output raster: C:\xxx\xxx\Documents\ArcGIS\Default.gdb\ifthe_ras's workspace is an invalid output workspace.

I'm stuck. I cannot understand what's wrong with my raster calculations.


texts: Do QGIS layer definition files (qlr) support joined layers?

texts:
I tried to save a qlr file for a shp layer joined to an xls file, but when I load it, the join definition isn't there anymore...

Join definition is present in the qlr file, so I assume it should work.

Has anybody already tried this ?


texts: ArcObjects: How to get a reference to layer in a java add in

texts:
I am trying to implement a java add-in version of this solution:
unique value symbolize layers from a lookup table.
At the same time I want to start coding in Java for the ArcGIS Desktop, so I am rewriting the sample Custom renderer, which is implemented as a property page in the layers -> properties -> Symbology -> Custom Renderer tab.
The above solution uses

IGeoFeatureLayer gfLayer = mxDoc.FocusMap.get_Layer(0) as IGeoFeatureLayer;    gfLayer.Renderer = (IFeatureRenderer)UvrItem.MakeUVR(@"Filepath");

Now, ArcMap surely knows the featurelayer, when editing its properties. But the Custom renderer is just handed a featureCursor. So I try to get a reference to GeoFeatureLayer
through the activeView like this:

    public boolean canRender(IFeatureClass fc, IDisplay display) throws IOException, AutomationException {

    try {
        HookHelper hookhelper = new HookHelper();
        IActiveView activeView = hookhelper.getActiveView();
        System.out.println("The Layer Name Is:   " + activeView.getFocusMap().getLayer(0).getName());
    } catch (Exception e) {
        System.out.println(e.getMessage() + e.getStackTrace().toString());
    }
    if (fc.getShapeType() == esriGeometryType.esriGeometryPolygon)
        return true;
    else
        return false;
}

But I get Error Message: Null.

If anyone more experienced than me knows the best way to do this in java,
I would very much appreciate some advice.

Best Regards, Marianne B. Wiese


texts: What is the fastest way to create a polygon of the part of PostGIS raster that is not filled with NoData values?

texts:
I have a raster table in PostGIS. For each row, I want to create a polygon that represents the part of the raster that is filled with data (not NODATA). I found out that ST_Polygon(rast) does the job. However, it works very very slowly.
Is there another solution for this problem that works faster?


texts:
I have a raster image which have to be aligned over the polygon.
Both the raster and polygon have the same coordinate system.

The figure shows the raster image which has to be fitted over the RED polygon.

https://drive.google.com/file/d/0B2RqG9tSAAIUc3J5ZDN6V01IcG8/view?usp=sharing

How can I do it guys?


texts: ERROR 000539: NameError: name 'prevX' is not defined

texts:
I'm trying to use the field calculator to input data into a field that I previously created, but it keeps saying that prevX is not defined in line 9. I'm very new to python and I don't really understand what's going on. This is my input into the code block

import math
def PtDistance(FID, Xval,Yval):
    global prevY
    global prevX
    if(FID == 1):
        prevX, prevY = Xval, Yval
        dist = 0
    else:
        dist = math.sqrt((prevX-Xval)**2+(prevY-Yval)**2)
        prevX, prevY = Xval, Yval
    return(dist)

and this is the error I keep getting

ERROR 000539: Traceback (most recent call last):
  File "< expression >", line 1, in <module>
  File "< string >", line 9, in PtDistance
NameError: name 'prevX' is not defined

Why is only line 9 the problem and not line 6 as well?


texts: Geoserver imagepyramid and Postgresql for granules

texts:
I'm running Geoserver 2.11.2 on Java 1.8.0_144 32 bits on Amazon Linux together with Tomcat 9. The plugins I have installed are imagemosaic, imagepyramid, GeowebCache (internal) and Geofence.

When using the imagemosaic plugin I am able to use a Postgresql DB in combination with JNDI connection pooling for storing the granules and time dimension via a datastore.properties file. The file looks like:

SPI=org.geotools.data.postgis.PostgisNGJNDIDataStoreFactory
schema=database_schema_name
Loose\ bbox=true
Estimated\ extends=false
validate\ connections=true
Connection\ timeout=10
preparedStatements=true
jndiReferenceName=java:comp/env/jdbc/postgis

Which work as expected. However when I add the same file in the imagepyramid datastructure, this setup fails. When I add this datastore.properties in just the lowest level (0) then that level gets configured in a database table, while the higher levels are configured using a shapefile. When I add a datastore.properties to all levels, geoserver returns the following error:

2017-09-28 07:39:40,780 ERROR [gce.imagemosaic] -
java.io.IOException
    at org.geotools.gce.imagemosaic.ImageMosaicConfigHandler.createGranuleCatalogFromDatastore(ImageMosaicConfigHandler.java:406)
    at org.geotools.gce.imagemosaic.ImageMosaicConfigHandler.createCatalog(ImageMosaicConfigHandler.java:329)
    at org.geotools.gce.imagemosaic.ImageMosaicConfigHandler.buildCatalog(ImageMosaicConfigHandler.java:962)
    at org.geotools.gce.imagemosaic.ImageMosaicConfigHandler.indexingPreamble(ImageMosaicConfigHandler.java:941)
    at org.geotools.gce.imagemosaic.ImageMosaicDirectoryWalker$MosaicDirectoryWalker.<init>(ImageMosaicDirectoryWalker.java:109)
    at org.geotools.gce.imagemosaic.ImageMosaicDirectoryWalker.run(ImageMosaicDirectoryWalker.java:197)
    at org.geotools.gce.imagemosaic.Utils.createMosaic(Utils.java:460)
    at org.geotools.gce.imagemosaic.Utils.checkSource(Utils.java:1336)
    at org.geotools.gce.imagemosaic.ImageMosaicReader.initReaderFromURL(ImageMosaicReader.java:532)
    at org.geotools.gce.imagemosaic.ImageMosaicReader.<init>(ImageMosaicReader.java:479)
    at org.geotools.gce.imagemosaic.ImageMosaicFormat.getReader(ImageMosaicFormat.java:504)
    at org.geotools.gce.imagepyramid.Utils.checkSource(Utils.java:200)
    at org.geotools.gce.imagepyramid.ImagePyramidReader.<init>(ImagePyramidReader.java:177)
    at org.geotools.gce.imagepyramid.ImagePyramidFormat.getReader(ImagePyramidFormat.java:249)
    at org.geotools.gce.imagepyramid.ImagePyramidFormat.getReader(ImagePyramidFormat.java:56)
    at org.geoserver.catalog.ResourcePool.getGridCoverageReader(ResourcePool.java:1511)
    at org.geoserver.catalog.ResourcePool.getGridCoverageReader(ResourcePool.java:1441)
    at org.geoserver.catalog.impl.CoverageStoreInfoImpl.getGridCoverageReader(CoverageStoreInfoImpl.java:59)
.....
    at org.apache.tomcat.util.threads.TaskThread$WrappingRunnable.run(TaskThread.java:61)
    at java.lang.Thread.run(Thread.java:748)
Caused by: java.lang.NullPointerException
    at java.lang.Class.forName0(Native Method)
    at java.lang.Class.forName(Class.java:264)
    at org.geotools.gce.imagemosaic.ImageMosaicConfigHandler.createGranuleCatalogFromDatastore(ImageMosaicConfigHandler.java:389)
    ... 160 more

Is it possible to configure the granules of all levels in the database? And if so, how would my datastore.properties file look like?
Like I said, it works fine for the lowest level (0) but it fails for all the levels above.

In case this is not possible, will it be sufficient to use the shapefiles as an index for the granules? The geoserver is used in a production environment and I am worried about performance issues when using the shapefiles as index.

Furthermore the imagepyramid automatically names the table after the level folder (0). This folder and table name are automatically created when configuring the datastore. Is there a way to custom name the folder and thus the table name? Otherwise I could end up with multiple '0' tables in the same schema.


texts: Garmin Fenix 6, map format for custom maps?

texts:
I was thinking about buying a new smartwatch and was looking at Garmin Fenix 6X PRO. Since I have the most detailed map available for Sweden (Fastighetskartan) in my handheld GPS (Garmin Alpha 100) I was hoping that there where some way of getting this map into the Fenix as well.
Does anyone know if there is some software (like QGIS, GDAL or other open source software) that could create the file type that Garmin Fenix can use?
And what file type it is that it requires?


texts: Changing page orientation using data driven pages

texts: Is it possible to change the page orientation of a map using data driven pages in ArcGIS 10.1? After setting up the DDP I have 80 portrait maps with all the map elements set up correctly. About 10 of the maps actually need to be in landscape, but with all the same elements, so when I am exporting them, I would like to be able to indicate that those 10 should be generated in landscape.


texts: Creating piechart of text column from attribute table in QGIS?

texts:
I would like to create a pie chart from a column of a attribute table, and I would like to see this pie chart above the map. The column contains names not numbers (I think fields aren't in a numeric format... on the other hand they're names!)
I work with marine species distribution data and I built a database of records of many many species...
In the specific, I have a column, called 'species', where there are many records (names) of some marine species. Several species may have many records, other species may have only a few numbers of records, so my objective is to graphically see the distribution of records among the species.
Please see the attachments:

img1 http://postimg.org/image/rn56c8l4z/
img2 http://postimg.org/image/e6918ynmj/


texts:
I am using Mapbox and I am trying to change the icon-image by clicking on it, for example I want to change this:

With this when I click:

map.on('load', function() {
    map.loadImage('./images/marker-removebg-preview.png',
    function(error, image) {
            if (error) throw error;
            map.addImage('marker', image);
        }
    );
});

var draw = new MapboxDraw({
    displayControlsDefault: false,
    controls: {
        polygon:true,
        trash: true,
    },
    userProperties: true,
      styles: [
          {
          "id": "gl-draw-polygon-stroke-active",//change le style des lignes de polygon
          "type": "line",
          "layout": {
            "line-cap": "round",
            "line-join": "round"
          },
          "paint": {
            "line-color": "#D20C0C",
            "line-dasharray": [0.2, 2],
            "line-width": 2
          }
        },

        {
          "id": "gl-draw-polygon-and-line-vertex-active",//change le style des points de polygon
          "type": "symbol",
          "filter": ["all", ["==", "meta", "vertex"], ["==", "$type", "Point"], ["!=", "mode", "static"]],
          "layout": {
            "icon-image": "marker",
            'icon-size': 0.06
            }
         }

       ]

    });var draw = new MapboxDraw({
    displayControlsDefault: false,
    controls: {
        polygon:true,
        trash: true,
    },
    userProperties: true,
      styles: [
          {
          "id": "gl-draw-polygon-stroke-active",//change le style des lignes de polygon
          "type": "line",
          "layout": {
            "line-cap": "round",
            "line-join": "round"
          },
          "paint": {
            "line-color": "#D20C0C",
            "line-dasharray": [0.2, 2],
            "line-width": 2
          }
        },

    {
      "id": "gl-draw-polygon-and-line-vertex-active",//change le style des points de polygon
      "type": "symbol",
      "filter": ["all", ["==", "meta", "vertex"], ["==", "$type", "Point"], ["!=", "mode", "static"]],
      "layout": {
        "icon-image": "marker",
        'icon-size': 0.06
        }
     }

   ]

});

I tried this method but it didn't work

map.on('click', 'gl-draw-polygon-and-line-vertex-active', function (e) {
    // load a smaller cat image
    map.loadImage('./images/marqueur_avec_camera.png', function(error, image) {
        if (error) throw error;
    map.addImage('marker-img', image);
    // change the icon-image property of the points layer to the smaller cat
    map.setLayoutProperty('gl-draw-polygon-and-line-vertex-active', 'icon-image', 'marker-img');

  });
});


texts:
I'm trying to have a thumbnail displayed and can get it working when using a specific picture path.

I have created a new column in the atribute table

I would like to have this be the image source in the code, but I can't seem to get it to work when putting in the name of the expression instead of a full path to a picture. Is there a way to tell the code to take the value of the table I am missing?

When I hover over a point I get:

With the following code:


texts: Checking if data has spatial definition?

texts:
I am trying to check if my data has a valid spatial reference and then check if its projected am running it an error saying can't execute:
My code:

import arcpy
import os, sys, traceback, shutil
import numpy as np

arcpy.env.overwriteOutput = True

try:
    outputfolder = arcpy.GetParameterAsText(0)
    fc = arcpy.GetParameterAsText(1)
    inputraster = arcpy.GetParameterAsText(2)

    #check to see if the arcinfo license is on
    if arcpy.CheckProduct("ArcInfo")=="Available":
        continue
    else:
        arcpy.AddError("An ArcInfo license is required to run this tool")
        sys.exit()

    #set the Progressor
    arcpy.SetProgessor("step","Checking input requirements....",0,9,1)

    #Process :create an out put folder
    arcpy.CreateFileGDB_management(outputfolder,"output.gdb")
    outputgdb = outputfolder + "\\output.gdb"+ os.step

    #create path to store raster data
    flooded = outputfolder +"floodedarea"

    #check to see if the inputs are in the same projection
    desccs_ras = arcpy.Describe(inputraster).SpatialReference
    desccs_fc = arcpy.Describe(fc).SpatialReference
    if desccs_ras.name,desccs_fc.name ="Unknown":
        arcpy.AddMessage("Your data does not have a defined spatial reference")
    else:
        if desccs_ras.projectionName == desccs_fc.projectionName:
            continue
        else:
            arcpy.AddError("The input Data should be in the same projection")
            sys.exit()
        del desccs_fc,desccs_ras

I don't know what is wrong .Please help.The error I get is:

SyntaxError: invalid syntax (zerowayRule.py, line 32)


texts: Is there a tool to draw GeoJSON points on a raster image?

texts:
I know there are tools like geojson.io which allow me to draw points on a real world map. But I work on a fictional map, so all I have is a big raster image. Later I want to show it in Leaflet using CRS.Simple.

So what I want is a tool that let's me put points on a raster image and gives me a GeoJSON with pixel coordinates instead of latlon.

My google-fu failed me, so I thought I'd ask here if anyone knows of a tool that can do this. Any ideas?


texts:
I installed QGIS through OSGeo4W, and I am trying to write a standalone PyQGIS script, however, I am unable to import the qgis, qgis.core, processing ect… modules, as shown in the picture, I get the ‘Import Error: No Module named qgis...

I sought help from the many examples that exist that require you to change your command line path to your python libraries, change your python path to your c:\OSGeo4W\apps\qgis\python ect…

Here is a list of the answers I tried, but, it seems myself and many others are still unable to fix this issue, as I could not get it to acknowledge the existence of those libraries:

Problem with import qgis.core when writing a stand-alone PyQGIS script:

Cannot run standalone QGIS script

Running Custom Applications:


texts: How to select within a distance using QGIS

texts:
With QGIS 2.12.3, I have two projected EPSG:2927 point shapefiles, A and B.  I would like to create a selection of A points that are within 100 meters of any B point.

This is easy in Arc, using the Select by Location tool with the search distance option.

I am aware of two Select by Location tools in QGIS (Vector -> Research Tools -> Select by Location and the Spatial Query plugin).  However, neither has a search distance option.

In summary, how to do a QGIS Select by Location with a search distance?


texts:
How to add waterways and railroads without turn restrictions into the previously created street network dataset  with turn restrictions, in order to create a multimodal network?

For a university project, we use TeleAtlas Multinet data. We used the following Street Data Processing Tool: http://www.arcgis.com/home/item.html?id=755f96fcde454ece8f790fecb3e031c7

Street Network is created with all kinds of turn restrictions, then now we need to add waterways and railroads into the routable network with global turns. However we couldn't create the multimodal network. Could you please give a suggestion if you have an idea about this topic.


texts: How to symbolize multiple layers with different EMF's in ArcObjects C#

texts:
I'm trying to symbolize multiple layers with an emf logo.  When debugging, it gets through the first layer then crashes.  I'm thinking it might have something to do with the variable "emf1" not changing.  I've tried a few things but can't get it to iterate properly.

    public void symbolize()
    {
        {
            IMxDocument pMxDoc;
            pMxDoc = (IMxDocument)ArcMap.Application.Document;

            IMap pMap;
            pMap = pMxDoc.FocusMap;

            IEnumLayer pLayers;
            pLayers = pMap.Layers;

            ILayer pLayer;
            pLayer = pLayers.Next();
            List<string> lyrName = new List<string>();

            while (pLayer != null)
            {
                if (pLayer.Name == "")
                {
                    break;
                }
                else
                {
                    lyrName.Add(pLayer.Name);
                }
            pLayer = pLayers.Next();
            }

            //QI
            IGeoFeatureLayer pGeoFLayer;
            pGeoFLayer = (IGeoFeatureLayer)pLayer;

            ISimpleMarkerSymbol pSym;
            pSym = new SimpleMarkerSymbol();
            pSym.Style = ESRI.ArcGIS.Display.esriSimpleMarkerStyle.esriSMSX;

            lyrName.ToArray();
            foreach (string name in lyrName)
            {
            IPictureMarkerSymbol emf1 = new PictureMarkerSymbolClass();
            emf1.CreateMarkerSymbolFromFile(esriIPictureType.esriIPictureEMF, "C:\\EMF\\" + name + "_Logo.emf");
            emf1.Angle = 0;
            emf1.Size = 24;
            emf1.XOffset = 0;
            emf1.YOffset = 0;

            ISimpleRenderer pRenderer;
            pRenderer = new SimpleRenderer();
            pRenderer.Symbol = (ISymbol)emf1;
            pGeoFLayer.Renderer = (IFeatureRenderer)pRenderer;
            }

            IActiveView pActiveView;
            pActiveView = (IActiveView)pMap;
            pActiveView.Refresh();
            pMxDoc.UpdateContents();
        }

    }

This is the error message that pops up, this is a pretty standard error message for me.


texts: why am I getting a horizontal solid strip when projecting large image?

texts:
I am using Mapnik to create tiles on a tile server.  I'm depicting sea level rise in North America using SRTM data, and it works fine in most of the continent.  However, when zoomed out in my browser (zoom levels 5 and lower), I get a peculiar horizontal stripe across the entire continent:

When I am zoomed farther in, this stripe does not appear, which makes me think it's not an issue with the data itself.  The stripe is from 51 degrees N to 52N.  I suspect it's some sort of projection issue, but I'm not sure how to troubleshoot.

Mapnik is being used with Tirex; it is reading from a .vrt file which in turn references ~2400 tifs.


texts: Merging multiple feature classes into one using Python

texts:
So I have feature classes from an ArcGIS database that I renamed using this Python code:

But now I have to merge all the feature classes that are output from this into one feature class named OSRS_ORN_NER. I'm not really sure how I should assign the inputs and outputs when I'm using the arcpy.Merge_management() tool. This is all I have right now:


texts:
Based on following information:

Size is 1200, 1400
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84", 6378137, 298.257223563, AUTHORITY["EPSG","7030"]],   (i parametri 2 e 3 sono il semiasse maggiore e lo schiacciamento)
            AUTHORITY["EPSG", "6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",42],
    PARAMETER["central_meridian",12.5],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre", 1, AUTHORITY["EPSG", "9001"]]
]
Origin = (-598884.437500000000000, 650790.250000000000000)
Pixel Size = (996.479919433593750,-999.677795410156250)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2018:04:01 00:10:06
  TIFFTAG_DOCUMENTNAME=/datamet/data/tmp/export/RADAR/36/RRN/VMI/VMI.geotif
  TIFFTAG_IMAGEDESCRIPTION=IDL TIFF file
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=IDL 8.3, Exelis Visual Information Solutions, Inc.
  TIFFTAG_XRESOLUTION=100
  TIFFTAG_YRESOLUTION=100
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -598884.438,  650790.250) (  4d32'13.20"E, 47d34'44.14"N)
Lower Left  ( -598884.438, -748758.664) (  5d56'18.16"E, 35d 4'36.70"N)
Upper Right (  596891.466,  650790.250) ( 20d26'12.34"E, 47d34'50.73"N)
Lower Right (  596891.466, -748758.664) ( 19d 2'23.69"E, 35d 4'40.94"N)
Center      (    -996.486,  -48984.207) ( 12d29'17.00"E, 41d33'32.31"N)
Band 1 Block=1200x1 Type=Float32, ColorInterp=Gray

I want to obtain (lon, lat) for each point of grid.

I'm trying with this Java code:

import ucar.unidata.geoloc.EarthEllipsoid;
import ucar.unidata.geoloc.LatLonPoint;
import ucar.unidata.geoloc.LatLonPointImpl;
import ucar.unidata.geoloc.ProjectionImpl;
import ucar.unidata.geoloc.ProjectionPointImpl;
import ucar.unidata.geoloc.projection.TransverseMercator;
import ucar.unidata.geoloc.projection.proj4.TransverseMercatorProjection;

public class TestGeoTiff {
    public static void main (String[] args) {
        //corners of my grid
        test(new double[] { -598884.438,  650790.250});
        test(new double[] { -598884.438, -748758.664});
        test(new double[] {  596891.466,  650790.250});
        test(new double[] {  596891.466, -748758.664});
    }

    public static void test (double[] XYproj) {

        ProjectionImpl proj = new TransverseMercator(42, 12.5, 1., 0., 0.);
        ProjectionPointImpl p   = new ProjectionPointImpl(XYproj[0], XYproj[1]);
        LatLonPointImpl lp      = (LatLonPointImpl) proj.projToLatLon(p);
        System.out.println(lp.getLongitude()+", "+lp.getLatitude());

        //or alternatively
        ProjectionImpl proj2 = new TransverseMercatorProjection(EarthEllipsoid.WGS84, 12.5, 42, 1., 0., 0.);
        ProjectionPointImpl p2  = new ProjectionPointImpl(XYproj[0], XYproj[1]);
        LatLonPointImpl lp2     = (LatLonPointImpl) proj2.projToLatLon(p2);
        System.out.println(lp2.getLongitude()+", "+lp2.getLatitude());
    }
}

The output is:

12.5, 90.0
-77.50000000000001, 1.229040709186457E-39
12.5, -90.0
-77.50000000000001, 9.003088929951159E-40
12.5, 90.0
102.50000000000001, 1.6804171192682352E-39
12.5, -90.0
102.50000000000001, 1.2309555453373566E-39

but it's not what I expected and it's very different from what it should be

Upper Left  ( -598884.438,  650790.250) --> (  4d32'13.20"E, 47d34'44.14"N)
Lower Left  ( -598884.438, -748758.664) --> (  5d56'18.16"E, 35d 4'36.70"N)
Upper Right (  596891.466,  650790.250) --> ( 20d26'12.34"E, 47d34'50.73"N)
Lower Right (  596891.466, -748758.664) --> ( 19d 2'23.69"E, 35d 4'40.94"N)

Is there a solution?


texts: Vertices editing / deletion question in QGIS (Blue vs. Red)

texts:
I am having an issue in QGIS.
I have two sets of lines in the same layer. The layer is in EDIT mode. Using the Vertex Tool, one line appears red (and I can use the tool), but the other line appears blue when highlighted with Vertex Tool, and I cannot move/edit the vertices (only delete).
Can someone explain why the second line does not appear red?


texts:
GetOpenFilename doesn't work in arcobjects VBA, so I have to use pGxDialog.

Standard pGxDialog doesn't let to open any file I want, a few extentions are available only (shp, dxf, xls...).

Maybe, I just don't know how to set a Filter.

I don't need to open extra extension files for reading/writing, I just need to return the PATH of the document to a variable, which will be used later on in the code.

I scripting in VBA, arcgis 10.2.2


texts:
I have following problem: I created a polygon shapefile. Now I want to label the polygons according to a certain column in the attribute table. No problem so far.
When I manually move the labels some labels disappear. The more labels I move the more labels disappear.
I tried many different things in the label settings as ticking "show all labels for this layer (including colliding labels)", "Label every part of multi-part features" or untick "Discouraging labels from covering features".
I have no idea how to solve that problem and would be happy with any ideas. I could also upload the data (but I don't know how to do it).


texts: How to select a line from a osm shape by osm_id

texts:
I'm using TileMill to make a map, but I can't seem to be able to format a specific item with a specific osm_id via mss.

I have a shapefile added as a layer #line:

Then, in the features list of #line, I can see the different osm_ids:

The problem is how can I format a specific line or element with a specific osm_id via the mss editor.

I tried this to make a very thick boundary lines for the city of Medellín:

#line[osm_id='1343264']{
  line-color:#555555;
  line-width:40;
}

and

#line[osm_id='1343264']{
  line-color:#555555;
  line-width:40;
}

As you can see, the element exists in OSM: http://www.openstreetmap.org/relation/1343264

AND it also exists in the #line shapefile I'm using. The osm_id in TileMill is the same, but with a minus sign in front.

How can I select this boundary line via mss to format it?


texts: Scatter plot in GEE with two series of data with different size

texts:
I am relatively new in GEE, and I am trying to plot what it should be a very simple scatter plot in another language such as R, but I am having a very hard time in GEE.
My goal is to display a scatter plot with two series with different size. Axes X and Y are representing Principal Components (PC) 1 and 2, respectively. The first serie represents values from flooded pixels and the second serie represents values from no flooded pixels. Therefore, I would like to display one serie: [PC1_Flooded,PC2_Flooded], and another serie: [PC1_NoFlooded,PC2_NoFlooded].
To distinguish flooded and no flooded pixels I selected a set of polygons on the GEE platform and set as FeatureCollection. Therefore, I have a geometry with several polygons for flooded, and another geometry with several polygons for no flooded.
The strategy I followed to get the scatter plot was to apply a mask with flooded geometry to PC1 and PC2 images, and a mask with no flooded geometry to PC1 and PC2 images. I used a mask since I did not want to distinguish among the different polygons in each geometry. The code I attach starts at the creation and application of the masks to the PC1 and PC2 images. Then, I merged PC1 and PC2 bands masked with flooded and no flooded geometries, resulting in an image with 4 bands: PC1_Flood, PC2_Flood, PC1_NoFlood and PC2_NoFlood, which contains all the data to represent the scatter plot.
//pcImage IS AN IMAGE WITH 4 PRINCIPAL COMPONENTS

//MASK CREATION FROM GEOMETRIES
var masked_flood = ee.Image.constant(1).clip(Flood.geometry()).mask()
var masked_noflood = ee.Image.constant(1).clip(No_Flood.geometry()).mask()

//APLICATION OF MASK TO PC IMAGE
var flood_masked_pc = pcImage.updateMask(masked_flood)
var noflood_masked_pc = pcImage.updateMask(masked_noflood)

//CREATE A NEW IMAGE WITH ONLY PC1 AND PC2 MASKED WITH FLOOD AND NO FLOOD GEOMETRIES
var pcImage_composite = flood_masked_pc.select(['pc1','pc2']).addBands(noflood_masked_pc.select(['pc1','pc2']))

//RENAME BANDS
pcImage_composite = pcImage_composite.rename(['pc1_Flood','pc2_Flood','pc1_NoFlood','pc2_NoFlood'])

//CONVERSION FROM IMAGE BANDS TO LIST OF VALUES
var pixelVals = pcImage_composite.reduceRegion({
  reducer: ee.Reducer.toList(),
  geometry: sesia,
  scale: 10,
  maxPixels: 1e13
})

//ERROR BECAUSE THEY HAVE A DIFFERENT SIZE
var yValues = pixelVals.toArray(['pc2_Flood','pc2_NoFlood'])

At this point is where I do not know how to proceed. I thought that somehow I could display the scatter plot with these 4 bands. I also read that I could convert these bands to a list of values, and then plot them.
I seek any help to solve this 'simple' problem. I followed this strategy, but I am sure that I could follow others such as retrieving values from geometries without applying masks.
Link to code: https://code.earthengine.google.com/fb928ca6124b5e1a9e9863c3763619ea


texts: Convertion of Geotiff from EU-DEM EPSG:3035 - ETRS89 / ETRS-LAEA to WGS84

texts:
I have downloaded elevation data from EU-DEM

My background is in software engineering and my GIS skills are limited.

How can I convert the aforementioned geotiff to WGS84 format geotiff with lon/lat values and altitude in meters.

I tried using QGIS "called warp command from the menu on the original file" however the converted lon/lat don't match with i.e. google earth as well as the altitude values definitely are not in meters.

Expanding my workflow as requested:

Ultimately I want to import and display elevation and other GIS data in an application made in Unity3D. My working environment is Windows.

The steps I followed are:

Convert specifically the "eu_dem_v11_E50N10.tif" (roughly the area convering Greece) file from QGIS with warp with the following settings:

Source CRS: empty

Target CRS: EPSG:4326 - WGS 84

No data value: default 0.0000

Output file res: default 0.0000

Resampling: Nearest neighbour

Advanced params: default

Output data: float32

Georeferenced extends: default

CRS of the target raster: default

Reprojected: c:\temp\eu_dem_wgs84.tif
Tried to read the file and display the elevation in a Terrain object within a Unity3D project using the library libtiff.net (Maybe here is the problem as I cannot read the heights correctly)

First issue is that, from my inexperience with QGIS, I did not find a way to find the altitude within QGIS of a lon/lat point.

Second issue is that using the gdalinfo on the generated file, the origin(lon/lat) is  wrong as in the original the origin is located north of Brindisi in Italy whereas in the converted the origin is located roughly between Bari and Taranto, resulting in reading wrong altitudes


texts: Identify most densely populated areas along streets with QGIS

texts:
On Google Satellite Imagery, I've created a layer with objects marking houses in a 1 km proximity to large streets in a rural setting using QGIS with the Openlayers Plugin. I would like to identify the top 10 most densly populated hotspots (most houses within a 1km circle) along app. 100 km accumulated street.

I am sure there is already much Information about such task on the web, it would already be very helpful to give me the right vocabulary for what I am trying to do to be able to do my own research.

Backround Information: My Task is to identify high demand spots for water in an underserved area in Western Africa. There is x number of water kiosks installed, connected by roads with an average distance of about 20km from one another. We would like to set up small water selling Points along these roads and therefore need to find suitable Locations which as many People as possible can reach by foot (i.e. 1 km max. distance).


texts: Which format to choose to convert to ShapeFiles?

texts:
Given the following format options:

all.json,
all.xml,
all.xls,
all.csv,
all.kml,
all.gjson

Which one is the best/easiest to convert the data into standard ERSI Shapefiles which contains the coordinates as well as the data attributes for each feature?

I tried importing KML and JSON via QGIS but the data attributes are squashed into a single column (see screenshot). Conversion via ogr2ogr was successful neither - it spits out a DBF file only. Conversion via csv2shp did not work at all - I could not figure out how to run the command.

I came up with the following all.vrt as suggested by afalciano.

<OGRVRTDataSource>
    <OGRVRTLayer name="all">
        <SrcDataSource>all.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>WGS84</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="laengengrad" y="breitengrad"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

All features are printed to the console running: ogrinfo -ro -al all.vrt


texts: Remove NoData from ArcGIS Raster Dataset

texts: I want to completely remove the NoData (pixel value 255) from a single band raster dataset. I don't want to make it transparent using symbology. Any ideas?


texts: Clipping raster image using mask layer is not clipping in correct way

texts:
I am using QGIS 3.10 when tried to clip raster with mask layer I got the same result of extent layer knowing that I activated option with extent boundary

Link for shapefile and raster


texts: osm2pgrouting how to configure mapconfig.xml file

texts:
I need to customize mapconfig.xml file.
I have to do bike routing so I start with mapconfig_for_bicycles.xml https://github.com/pgRouting/osm2pgrouting/blob/master/mapconfig_for_bicycles.xml
This configuration ignore, for examples, ways like
https://www.openstreetmap.org/way/409033340
with tag name "bicycle" and tag value "designated" (but also tag name "highway" and tag value "path").
So I edit the xml file adding this rows:
....
<tag_name name="bicycle" id="5">
<tag_value name="designated" priority="1.6" id="501" />
</tag_name>

<tag_name name="highway" id="6">
<tag_value name="cycleway" priority="1.0" id="601" />
<tag_value name="path" priority="1.1" id="602" />
<tag_value name="pedestrian" priority="1.1" id="603" />
...

With this configuration after import, tag_id is still  602, not 501 as expected.
To have expected result (tag_id = 501) I need to delete the row:
<tag_value name="path" priority="1.1" id="602" />

But in this way I'm not sure I will import all ways with tag name "bicycle" and tag value "designated" and I cannot have a tag_id for all ways with tag name "highway" and tag value "path".
What is the best way to configure xml file to include the 2 info together?


texts: What's the unit of measurements produced by "sum line lengths"?

texts:
My goal is to find the total length of roads per county for a particular state. You would think there would be a simple government spreadsheet with this information in it....

I cannot figure out how to find the summed lengths of a polyline layer(roads) in each county (polygon) in QGIS. I have tried using Vector-> Analysis -> Sum line lengths but it just reproduces the county shapefile with a length column, of which I have no idea what unit it is in.


texts: Sort Custom Tileset Features in MapBox Studio or GL

texts:
I have a large number (hundreds of thousands) of point features uploaded to my account in MapBox studio as a dataset, and have generated a tileset based on that dataset.  This tileset is used in a Style as a layer source for a circle layer.  I am then using this style as the basis for a slippy map implemented in MapBox GL JS
How do I use one of the properties of the features in my datasource to determine the drawing order of the features in mapbox gl?  I can't see any options for that in the style editor. I've tried setting the circle-sort-order layout property
map.setLayoutProperty('<layername>','circle-sort-key',["-",["get","<propertytosorton>"]])
but that did not seem to change anything.


texts: QGIS GRASS vector simplification creates duplicate data rows

texts:
When simplifying polygons in QGIS (2.2.0) with the GRASS v.generalize tool the resulting layer duplicates some data rows many times.  A shapefile of 657 zones/data rows:

.. is simplified to one with correct-looking geometries but 12,513 data rows:

Many of these don't seem to have an associated polygon, at least not a visible one.  Not sure why this happens, but what would be the easiest way to (preferably) prevent it or (alternatively) clean the data table up?

The problem is replicable using this small shapefile with CRS set to EPSE 27700 - 33 data rows become 60 after generalising by the 'douglas' method with a 1000 tolerance.

Many thanks in advance.

Edit.

I suspect this has something to do with the clipped nature of the shapefile which follows coastlines and river edges (which do not all simplify as you might predict in line with the tolerance).


texts: GEE Process for deciding what pixels lie within a geometry

texts:
I am using image and collection reducers to extract mean values of several assets (for example, prism precipitation and pml2 ET) from within USGS watershed boundaries which I load in as geometries. I am generally extracting data at 1000m scale.
I am having no problem with this workflow, however I am wondering what GEE does to determine which pixels qualify as within the geometry? (Is it based on the centroid or edge or another calculation?)