texts: How to represent two fields of the same states or district layer using different color borders and fills in QGIS?

texts:
If I had a set of data in QGIS that was gradually shading a state or district the exact way I wanted -- but I also need to portray more information about that state or district with additional colors, is that possible to do using colored borders? If so, can somebody walk me through that? And if not, has anybody else thought of another way to represent multiple things about one location?

I guess the ideal would be something like shading the state lame yellow to bright yellow pending on how high their statistic was -- and if the state was Democratic-controlled, put a blue border around it, and if it was Republican-controlled, put a red border around it. Something like that.

Any suggestions? Thanks!

--

*Edit: I think I've potentially figured it out by adding another layer and creating my own border-only style and classify one as Democrat (blue) and one as Republican (red). If the polygon fills are 100% transparent and this layer was on top of the other one, I'd imagine that this should work, right? But every time I get to hitting "ok" to where the map should change to reflect the new borders, it always loads forever to the point where I have to 'Force Quit.'

Any thoughts on that? Thanks again!


texts:
I'm new to using GIS data.

Using QGIS I'm trying to export a DXF file from a .shp file from the NYC OpenData site.

I can open the file, pan, zoom, select the area I want to export, etc, but when I go to Project>DXF Export the export fails every time.

There are many different settings within the DXF Export dialog box, but I'm not sure what they mean or what should be selected.

The map data I'm trying to export are very simple NYC building footprints.

I'm on a Mac running OS10.10.

Ideally I want to export a DXF file that I can import into SketchUp.

Can anyone help with this?


texts: Geoserver polygon labeling Arabic numbers using SLD

texts:
I've used SLD to label polygon features with their parcel Ids. The problem is for Ids include '.' and '-' marks, the ordering of character is messed up. For example 39.1 changes to 139. I've checked them in Pgsql and ArcMap which show them correctly.

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:se="http://www.opengis.net/se" xmlns:xlink="http://www.w3.org/1999/xlink"
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.1.0"
  xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.1.0/StyledLayerDescriptor.xsd">
  <NamedLayer>
    <se:Name>parcels</se:Name>
    <UserStyle>
      <se:Name>parcels</se:Name>

      <se:FeatureTypeStyle>
        <se:Rule>
          <se:Name>Single symbol</se:Name>
          <se:MaxScaleDenominator>2000</se:MaxScaleDenominator>
          <se:PolygonSymbolizer>
             <se:Stroke>
              <se:SvgParameter name="stroke">#7e7e7e</se:SvgParameter>
              <se:SvgParameter name="stroke-width">0.7</se:SvgParameter>
              <se:SvgParameter name="stroke-linejoin">bevel</se:SvgParameter>
            </se:Stroke>
          </se:PolygonSymbolizer>
          <se:TextSymbolizer>
            <se:Geometry>
              <ogc:Function name="centroid">
                <ogc:PropertyName>geom</ogc:PropertyName>
              </ogc:Function>
            </se:Geometry>
            <se:Label>
              <ogc:PropertyName>id_fa</ogc:PropertyName>
            </se:Label>
            <se:Font>
              <se:SvgParameter name="font-family">Iranian Sans</se:SvgParameter>
              <se:SvgParameter name="font-size">11</se:SvgParameter>
              <se:SvgParameter name="font-weight">normal</se:SvgParameter>
            </se:Font>
            <se:LabelPlacement>
              <se:PointPlacement>
                <se:AnchorPoint>
                  <se:AnchorPointX>0.5</se:AnchorPointX>
                  <se:AnchorPointY>0.0</se:AnchorPointY>
                </se:AnchorPoint>
                <se:Displacement>
                  <se:DisplacementX>0</se:DisplacementX>
                  <se:DisplacementY>5</se:DisplacementY>
                </se:Displacement>
              </se:PointPlacement>
            </se:LabelPlacement>
            <se:Fill>
              <se:SvgParameter name="fill">#523735</se:SvgParameter>
            </se:Fill>
            <se:VendorOption name="autoWrap">150</se:VendorOption>
            <se:VendorOption name="maxDisplacement">50</se:VendorOption>
            <se:VendorOption name="group">yes</se:VendorOption>
          </se:TextSymbolizer>
        </se:Rule>
      </se:FeatureTypeStyle>

    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

I should mention that I've translated English numbers to Arabic format using python and inserted them in a new field. Is it something wrong with that process?


texts: Export Mapbox-Studio project to Mapnik XML

texts:
Tilemill had a feature to export the project to Mapnik XML (see here). Is there an equivalent functionality in Mapbox Studio?

Some background: I'm trying to use a style like OSM Bright, edit it a little and the use in on my own tileserver (running mapnik).


texts: Replace old invalid polygons with new fixed ones

texts:
I have a layer of 2053 polygons. Some of these have errors.
I corrected these errors with Processing LWGEOM Provider.
The outputfile have the corrected polygons only.

How do I merge the corrected ones with the old and in the same time deleting the old polygons with errors?

I've attached 2 pictures.
The Green layer is the corrected polygons
The pink layer is the original layer with some errors


texts: Looking for an algorithm to quickly find the closest point on a road network (or other graph)

texts:
I am looking for a fast algorithm that will find the closest point on a line segment in a graph (such as a road network) to my test point.

I'll try to be more specific. I have a road network graph. This may have somewhere in the region of 30 million nodes and links/lines (my US dataset is of this size). I already have a fast routing algorithm up and running, but it works from node to node. The real world does not work with nodes, but with coordinates. I could search for the nearest node to each start/end coordinate, and in the past I have coded a very fast coordinate matching algorithm (a few million nodes, and managing >~1000 test points per second).

Commercial routing algorithms do not work like that: They  direct you to the nearest road segment. In other words I need to match a coordinate to the nearest line (graph link). As well as finding the closest line, I need a fraction that tells me how far along the line that closest point is. Then I can route from the closest node and add the fractional line in to the distance/time sum. Or correct accordingly if the route happens to go down the line.

I think it is reasonable to have a limit on the maximum distance to search. E.g. "The closest line that is within 500 metres". And this should help prune the search space.

Coding up an algorithm that searches all of these lines is pretty simple, but it would be slow due to the large number of lines and the geometry calculations. I could limit the search to be limited to those lines within the distance limit. And this could be further optimized by sorting the lines by longitude (say). Of course most lines have a fine longitude "width" which complicates that sort&search.

I still think this will be slow. How do commercial packages do this? Most will find and display a route within half a second or so, and most of that is taken up with the main routing algorithm, and constructing the result (deriving turn directions, etc).

There are a number of similar questions on this site (eg. How to find the nearest point projected on the road network? ) but these refer to specific packages such as PostGIS or ArcGIS. I am looking for an actual algorithm to code up myself. I suspect my main issue is known the words/phrases to search for. Once I've found an algorithm name, then a bit of Googling/etc will quickly find all-sorts of algorithms and papers...


texts: Can I make Points-in-polygon output to the attribute table of polygon layer?

texts:
I have three layers of points, corresponding to buildings of different types, on top of polygons corresponding to output areas.

Running the points-in-polygon function seems to make a new polygon object with the pointcount in its attribute table

I'd like to instead have the pointcounts as new columns in the output area polygon layer, so that I can colour them based on dominant building type.


texts: Developing GUI in Python for ArcGIS geoprocessing using PyQT/Tkinter/wxPython?

texts:
I want to develop a GUI in Python for ArcGIS geoprocessing.

Can I use PyQT for GUI programming in ArcGIS?

I have also considered using Tkinter and wxPython.


texts: Comparing identical polylines direction in ArcMap

texts:
I have a file with polylines digitized in several directions (clockwise and counter-clockwise) and I need them to be all digitized in only one direction.
This I have solved with creating coverages and then back to lines again. Then I have dissolved them by Z value to get rid of short line pieces or duplicates that somehow appear.

Now I would like to find out actually how many of those new lines have been "fixed" - direction changed from false to correct.
Feature compare works on simple and exact duplicates (geometry change), but when the line has been merged for example, then it immediately gives geometry error or as is with my previous solution - the endpoint isn't the same any more and geometry is different.

Is there a solution?

Polylines are contours. More precisely depression contours. They are represented by "ticks" showing the direction of descent. The system representing them is locked, so the depressions have to be made according to this, but errors happen... a lot.
Left is digitized counter-clockwise and right is digitized clockwise. Left is correct.

I haven't found any tool that could find which closed contour is digitized wrong or correct, only thing that has helped is the conversions.
But now I need to start finding how many of these errors occur and I thought that if I can't find the errors themselves, I can find how many objects have changed.


texts: I have seen am ability in SQL Server 2008 Spatial Tools here that saves result of query as shapefile..I am intrested in how to do this task using C#..Is there any way to do this?I found ways to convert WKT to WKB and so on but I could not figur out how to convert them to .shp.I found a tool called ogr2ogr but I dont know how to use it in my C# code


texts: Computing distance from points to coastline

texts: I have a shapefile with the nodes of a grid I have applied, for the purposes of interpolation, and I want to compute the euclidean distance that each grid node (it is a point vector) has to the coastline. I don not want to calculate the nearest neighbor for each point of my grid. I want to have in a table the euclidean distance that each point (the points are placed on the nodes of a regular grid) has from the coast, so that I can use this information as an auxiliary variable in my interpolation. I have found several tools in ArcMap that compute distances, but nothing that computes all the distances from one point feature to a line feature.


texts: bringToFront() not working on tooltips in Leaflet

texts:
I have a Leaflet map that includes polygons and data from a GeoJSON file. I have permanent tooltips for each polygon that include data for the user to see. But there are so many tooltips that sometimes they overlap, and then the reader can't read the overlapped tooltips.

A solution to this problem, I thought, would be that if the user clicked a tooltip, it would be brought to the front. Unfortunately, I haven't been able to make this idea work. Does anyone know how to make this happen please?

Here is my code:

//define the map
var map = L.map('map').setView([40.350641, -79.86773], 15);

var OpenStreetMap_France = L.tileLayer('https://{s}.tile.openstreetmap.fr/osmfr/{z}/{x}/{y}.png', {
    maxZoom: 22,
    ext: 'png',
    attribution: '&copy; Openstreetmap France | &copy; OpenStreetMap contributors'
}).addTo(map);

//grab data about real estate parcels from json file and put in tooltips
var parcels = new L.GeoJSON.AJAX("data/parcels.geojson",{ //add parcel shapes and data
    onEachFeature: function (feature, layer) {
        layer.bindTooltip(/*feature.properties.PIN + " " + */"
" + feature.properties.PROPERTYHOUSENUM__asmt + " " + feature.properties.PROPERTYADDRESS__asmt + "
" + "
Use in 1922
" + feature.properties.USE_1922 + "
" + "
Current Use
" + feature.properties.CURRENT_USE, { 'permanent': true, 'interaction': true });

        //when you click a tooltip, bring it to the front
        layer.on('click',function() {
           this.bringToFront(); //doesn't work
           this.setZIndexOffset(999); //doesn't work
        });

    }
}).addTo(map);


texts: Looking for Parish Boundaries Datasets for QLD Austraila

texts: Every so often, I am asked for a parish boundary dataset. I can never find them as I am of the belief that they are no longer used. Has any one had any experience with parish boundaries in Australia (related to the state of Queensland) that could shed some light on the subject. I do know that the current Cadastral (DCDB) does contain the parish attributes (see answer below), but  I guess what I am really looking for is a current parish boundaries dataset (if it exist) that I could use to check against Cadastre (DCDB).


texts:
I have a geocoded raster image that I want to use as the background to a vector map with D3.js. The screenshot shows the layers in QGIS (grey is raster). Here's the markup for the map:

<div class="container"> <-- has raster image as background-image -->
    <svg></svg>
</div>

Also, see this map from the New York Times for an example of what I'm aiming for. My question is: How do I make the vector layers and the background image align in the browser?

I think this is the relevant code sections from the NYT, but I'm not sure how to apply it to my own project:

var background = $('.the-background');
var basemapAspectRatio = 0.6244;
var mapCenter = [42.45,32.42];
var mapScale = 2.288;
var mapRatio = 7000 / 4371;
var mapHeight = background.height();
var trueWidth = mapHeight * mapRatio;
var mapWidth = background.width();

projection = d3.geo.mercator()
  .center(mapCenter)
  .scale(trueWidth * mapScale)
  .translate([mapWidth / 2, mapHeight / 2]);


texts: How to modify the edit widget in ESRI flex viewer 3.0?

texts:
I am running ArcGIS Server 10.0 (sp2) on Windows server 2008 64-bit with IIS 7.  I am using the ESRI flex viewer version 3.0.

I have created a feature layer and published it as a rest service and included a template and use the edit widget with this feature service to allow specific users to modify a single attribute on road features on the feature layer.  The Edit widget, as is, will allow you to select multiple features (in this case, roads), but you cannot update all selected features at once, you have to do it one at
a time.  I want to add a button to the window that appears (the attributes window) that will allow you to apply the change to all selected features.

I have flex builder 4.6 and have created some custom widgets, but this change is a little beyond my skillset currently.  I've looked in the editwidget.mxml and editwidgetattributeinspectorskin.mxml, but don't see where actual edits are made in either of these files.  I'm thinking the actual codeing for applying edits to the selected features is in a different section of the source, but don't know where to look.

Any help is greatly appreciated.


texts: ST_Clip tiled raster using geometry results in out of geometry pixels

texts:
I'm trying to run an ST_Clip for a simple AOI geometry.

I'm using this WGS84 raster, imported into PostGIS like so:

raster2pgsql -t 100x100 -C -M -F -I -s 4326 Hansen_GFC-2016-v1.4_lossyear_20S_070W.tif hansen_temp | psql

I'm trying to clip this raster to my AOI (gist here).

I use ogr2ogr to insert this AOI into postgres:

ogr2ogr -f Postgresql PG:"dbname=charlie" aoi.geojson

I run the following SQL:

CREATE TABLE ras_clip AS
SELECT *
FROM (SELECT ST_Union(ST_Clip(rast, wkb_geometry ))
            FROM hansen_temp, aoi
            WHERE ST_Intersects(rast, wkb_geometry)) foo;

This results in a raster that seems to follow my AOI - perhaps following the 100x100 tiles I set when importing the data using raster2pgsql. Here's a screenshot taken from PostGIS:

Is this the desired results of ST_Clip? Is there any way to constrain the output to only pixels that actually intersect the geometry of interest?

PostGIS version info:
POSTGIS="2.4.4 r16526" PGSQL="100" GEOS="3.6.2-CAPI-1.10.2 4d2925d6" PROJ="Rel. 5.0.1, April 1st, 2018" GDAL="GDAL 2.3.0, released 2018/05/04" LIBXML="2.9.7" LIBJSON="0.13.1" RASTER


texts: Using Field Calculator with Python Parser and if/elif statements in ArcGIS for Desktop?

texts:
I have tried all of below snippets to use Python Parser in Field Calculator and update the values of a field called type based on a filed called MamerMN but in all of them I am getting Syntax error in Geoprocessing result window!

if !MamerMN! <= 0.151560:
    return 1

and

if (!MamerMN! <= 0.151560):
    return 1

and

if (MamerMN <= 0.151560):
    return 1

and

def(MamerMN)
if MamerMN <= 0.151560:
    return 1

and

def(MamerMN)
if (MamerMN <= 0.151560):
    return 1

Can you please let me know what I am doing wrong?

Update

def foo(MamerMN):
  if MamerMN <= 0.151560:
     return 1


texts: Repairing disconnected layers with selection sets using ArcObjects?

texts:
Has anyone had any success of repairing disconnected layers that have selection sets that have their data sources broken in ArcGIS?  Is there a way to detect this through code and then re-create that selection set?

I tried:

IFeatureLayerDefinition.DefinitionSelectionSet.Count

and

IFeatureSelection.SelectionSet.Count

of which just give an "Object reference not set to an instance of an object" error.

We are on version 10.0 and code in VB.NET or C#.


texts: Using LIKE clause in update cursor gives syntax error?

texts:
I am trying to write a second update cursor within one script.  This requires a LIKE condition to look up a value within the same layer to identify rows to update the field within that layer (first update cursor looks up a second table and is working successfully).  When I try to run the query it is throwing a syntax error and highlighting LIKE.

See my initial question about Writing python script to select attributes not in growing table with over 4000 values?.

Here is the portion of script in causing the issue:

with arcpy.da.UpdateCursor("in_memory/headphotos",[field2,'Located']) as UCur1:
    for URow1 in UCur1:
        if URow1[0].lower() LIKE '%reserve%':
            URow1[1] = 1

Here is my full script:

import os, sys, arcpy

arcpy.env.workspace = (directory location)
jpgList = arcpy.ListRasters('*','JPG')
photos = ([s.lower().strip('.jpg') for s in jpgList ])

arcpy.CopyFeatures_management(feature class location,r"in_memory/headphotos")

arcpy.AddField_management(r"in_memory/headphotos",'Located','SHORT')

with arcpy.da.UpdateCursor("in_memory/headphotos",['ce_id','Located']) as UCur:
    for URow in UCur:
        if URow[0].lower() in photos:
            URow[1] = 1
        else:
            URow[1] = 0
        UCur.updateRow(URow)

with arcpy.da.UpdateCursor("in_memory/headphotos",[field2,'Located']) as UCur1:
    for URow1 in UCur1:
        if URow1[0].lower() LIKE '%reserve%':
            URow1[1] = 1

arcpy.MakeFeatureLayer_management("in_memory/headphotos","tmp_layer",'Located = 0')
arcpy.CopyFeatures_management("tmp_layer",new layer)


texts: Calculating sink volume in R via rgrass7

texts:
I want to calculate the volume of sinks in a DEM raster in R, since I have to automatize the workflow later. Following the results of my previous question (Creating depressionless DEMs with rgrass7) I have my rasters clipped to the extent of the sinks. I did that by doing:

sinks <- regular_raster.tif - depressionless_DEM.tif
sinks[sinks == 0] <- NA

Now I'm looking for a way to calculate the volume of all the sinks in this raster. In ArcGIS or QGIS I would region group the raster, then calculate a contour file and a TIN for every sink and calculate its volume.
Within GRASS, all I find are ways to calculate a DEM from a contour file (like the v.surf and r.surf packages), but not the other way around. What am I missing? Or is the workflow in GRASS different?


texts: Search radius for some intermediate distance

texts:
ArcMap on my laptop runs very slowly. So it can take a long time for a near table to run.

Is it possible to use a search radius for some intermediate distance?

For example, the search radius only lets me specify up to 2000 metres, or up to 4000 metres.

My question is, is it possible to calculate distances in the range 2000-4000 metres?


texts: constrain random points with a polygon in R using sp library

texts:
I have a map of London's waterways, which I've grided into a raster dataset to enable me to measure proximity to water from randomly generated points.

What I'd like to do please is to constrain the randomly generated points to a polygon (which represents the sampling area) within  the larger grid.

This is the code I've used:

randLocs <- spsample(waterGrid, n=200, type='random')

waterGrid is the overall grid enabling me to measure distances, and I have a polygon which I've loaded as a shapefile.


texts: Checking for exact value in QGIS raster calculator

texts:
I would like to single out values which are equal to a certain floating point value using QGIS raster calculator. I tried to run the following:

"Converted@1"==-3.40282347*10^38

But I got the error message:

  Incorrect parameter value for LAYERS

Why does this happen? How can I check for an exact value in QGIS?


texts: Publish SQL Views from Postgres/Postgis to Geoserver

texts:
I've trying and reading how to publish a Postgres/Postgis View in Geoserver. I know that I can create an SQL View in Geoserver, but I'd like to define the View in Postgres, not in Geoserver.

Does anyone has an example of how it's done?

So far I've tried this, but I can't see the content of the View, not even in Postgres.

CREATE OR REPLACE VIEW camino_view AS
SELECT c.id, c.nombre, t.id AS tipo_id, t.nombre AS tipo, l.gid AS geo_id, st_transform(l.the_geom, (-1)) AS the_geom
FROM camino c, camino_tipo t, camino_layer l, camino_link e
WHERE c.tipo_id = t.id AND e.camino_id = c.id AND e.camino_geo_id = l.gid;

CREATE TABLE camino_layer
(
gid serial NOT NULL,
"Tipo" character varying(50),
"Nombre" character varying(50),
the_geom geometry,
CONSTRAINT camino_layer_pkey PRIMARY KEY (gid),
CONSTRAINT enforce_dims_the_geom CHECK (st_ndims(the_geom) = 2),
CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'LINESTRING'::text OR the_geom IS NULL),
CONSTRAINT enforce_srid_the_geom CHECK (st_srid(the_geom) = (-1))
)
WITH (
OIDS=FALSE
);

CREATE TABLE camino_link
(
 camino_id integer NOT NULL,
 camino_geo_id integer NOT NULL,
 CONSTRAINT camino_link_pkey PRIMARY KEY (camino_id, camino_geo_id),
CONSTRAINT camino_link_camino_fkey FOREIGN KEY (camino_id)
  REFERENCES camino (id) MATCH SIMPLE
  ON UPDATE NO ACTION ON DELETE NO ACTION,
CONSTRAINT camino_link_camino_layer_fkey FOREIGN KEY (camino_geo_id)
  REFERENCES camino_layer (gid) MATCH SIMPLE
  ON UPDATE NO ACTION ON DELETE NO ACTION
)
WITH (
OIDS=FALSE
);

CREATE TABLE camino_tipo
(
id serial NOT NULL,
nombre character varying(48) NOT NULL,
CONSTRAINT camino_tipo_pkey PRIMARY KEY (id)
)

CREATE TABLE camino
(
id serial NOT NULL,
nombre character varying(40) NOT NULL,
tipo_id integer NOT NULL,
CONSTRAINT camino_pkey PRIMARY KEY (id),
CONSTRAINT camino_tipo_fkey FOREIGN KEY (tipo_id)
  REFERENCES camino_tipo (id) MATCH SIMPLE
  ON UPDATE NO ACTION ON DELETE NO ACTION,
CONSTRAINT camino_name_uq UNIQUE (nombre)
)
WITH (
OIDS=FALSE
);

Any help will be appreciated.
Thanks,


texts: QGIS crashes when I try to check for double or check for crossing lines

texts:
I have a LineString layer that I want to check for double geometries and/or intersecting lines. I've tried two different tools (Check Geometry (Geometryprüfung) and Check Topology (Topologyprüfung)) and they both crash with my layer and I have no idea why.
Crash resport:
Crash ID: 12b1bee7ec347015efdb44f06b72d8dacdfe3999

Stack Trace

Index::~Index :
Index::~Index :
Index::~Index :
QgsSpatialIndex::addFeature qgsspatialindex.cpp:453
QgsSpatialIndex::addFeature qgsspatialindex.cpp:410
QgsFeaturePool::insertFeature qgsfeaturepool.cpp:126
QgsVectorDataProviderFeaturePool::QgsVectorDataProviderFeaturePool qgsvectordataproviderfeaturepool.cpp:41
QgsGeometryCheckerSetupTab::runChecks qgsgeometrycheckersetuptab.cpp:445
QMetaObject::activate :
QAbstractButton::clicked :
QAbstractButton::click :
QAbstractButton::mouseReleaseEvent :
QWidget::event :
QApplicationPrivate::notify_helper :
QApplication::notify :
QgsApplication::notify qgsapplication.cpp:458
QCoreApplication::notifyInternal2 :
QApplicationPrivate::sendMouseEvent :
QSizePolicy::QSizePolicy :
QSizePolicy::QSizePolicy :
QApplicationPrivate::notify_helper :
QApplication::notify :
QgsApplication::notify qgsapplication.cpp:458
QCoreApplication::notifyInternal2 :
QGuiApplicationPrivate::processMouseEvent :
QWindowSystemInterface::sendWindowSystemEvents :
QEventDispatcherWin32::processEvents :
UserCallWinProcCheckWow :
DispatchMessageWorker :
QEventDispatcherWin32::processEvents :
qt_plugin_query_metadata :
QEventLoop::exec :
QCoreApplication::exec :
main main.cpp:1637
BaseThreadInitThunk :
RtlUserThreadStart :

QGIS Info
QGIS Version: 3.16.5-Hannover
QGIS code revision: 58ba7c1ed6
Compiled against Qt: 5.11.2
Running against Qt: 5.11.2
Compiled against GDAL: 3.1.4
Running against GDAL: 3.1.4

System Info
CPU Type: x86_64
Kernel Type: winnt
Kernel Version: 10.0.18363

I've also tried on two different machines, exported the layer and projected it but nothing helped.
Does anybody have a solution?


texts:
If I upload my KML file, which contains polygons and pins, to Google Maps, then most of the polygons are missing (see below).

To try and address this, I pasted the kml url to the Google Maps search bar and all of the polygons are now shown on the map but all the pins are missing (see below).

The KML file can be downloaded from here.

What should I modify in the KML source to have all polygons and pins in the same map?


texts: Incidence angle of geographic cordination

texts:
I am processing Sentinel 1 data with ENVI SARscape (can process also with gamma if it helps me) for earthquake inversion. Due to the fact that I want to use ascending and descending acquisition, I am working with the earth coordination system and not the satellite one. Therefore I need to know the incidence angle of the ascending and descending track for specific geographic coordination points.

So I think my question is: How do I convert the geographic coordination to the satellite coordination system (for ascending or descending)?

After knowing the coordination of a point in the satellite coordination system calculating the incidence angle is clear.


texts: Enhancing/smoothing Landsat 8 image using ArcGIS Online?

texts:
I understand that Landsat 8 is 30m resolution, but is there any technique that can be used to smooth/enhance the NDVI or Agriculture Landsat layers in an ArcGIS Online (AGOL) web map?

I have adjusted the image quality to 'best' but can I achieve better results than this for "Agriculture" or "Agriculture DRA" processing template;


texts: Altimetric image with blank space when merging

texts:
I have downloaded some tiff images that has altimetric data from the Brazilian SRTM website: where i am using raster layers to apply merging. However, i have some problems with my merged images, a portion of my map is blank:

These areas correspond to regions: SC-20-zc, SD-20-XA, SC-20-Z-D, SD-20-X-B, SD-20-X-C, SD-20-X-D.
However when i plot them inside the state/province limits just using the raster added layers, they show up just fine, why is this happening?
I've tried visualizing the raster layers of that area that has bright color/blank space and this is the output:

This is my hypsometric map:

The light blue area is just missing, is this something wrong with how i performed the merging?
This is the actual map where apparently it's the yellow zone that represents the texture portion is blank in my map:

I'd like to know, what did i do wrong and how do i fix it?


texts: Not able to get lat in ol.proj.transform

texts:
I'm not sure what's going on with this code, but I don't see anything wrong with it.

My alert is able to return the lon and lat exactly as I need them. However, when I pass it through ol.proj.transform and ol.geom.Point, lat is never returned.

This is my HTML:

Virginia, United States

This is my javascript:

$(document.body).on('click', '.query-result', function(){
    var lon = $(this).attr('data-long');
    var lat = $(this).attr('data-lat');
    alert(lon);
    alert(lat);
    var feature = new ol.Feature(
      geometry: new ol.geom.Point(ol.proj.transform([lon, lat], 'EPSG:4326', 'EPSG:3857'))
    );
    feature.setStyle(iconStyle);
    vectorSource.addFeature(feature);
    map.addLayer(vectorLayer);
});

I'm specifically plotting a point on map when someone clicks .query-result. It works when I hardcode the coords, but not when passing them through my code.

I don't understand why.


texts: Display a messageBox with an image in ArcMap

texts:
I wanted to add "ABOUT" button to a Pythonaddins toolbar
the button should display a window that contains both of a text defines the authors and an image about the developer organisation.
I followed an example about using Tkinter, as the following adress
http://gisstudycenter.blogspot.com.eg/2014/12/how-to-create-graphical-user-interface.html

the script works as shown

but immediately the arc crashes

The code I used is

from Tkinter import *
import datetime
from array import *

root = Tk()
logo = PhotoImage(file="D:\My Papers\Paper 19 - StandAlone Gis Application\Draft\William_Shakespeare.gif")
w2 = Label(root, image=logo).grid(row=0,column=1)
explanation = """Team Geo Zoner"""
w2 = Label(root, compound = CENTER, fg="Red", bg="black", font="Times 55 bold", text=explanation)
w2.grid(row=0,column=0)

p=str(datetime.datetime.now().isoformat())
texter="Your Heading Goes Here"+str(p);
w3 = Label(root, compound = CENTER, fg="Green", font="Times 20 bold", text=texter)
w3.grid(row=1,column=0)

Parameters = ['Number of Inputs','Number of Ouputs','Number of Processing Parameters','Number of Intermediate data']

r = 2
i=1
for c in Parameters:
    Label(text=c, relief=RIDGE,width=25).grid(row=r,column=0)
    Entry(text="Value", relief=SUNKEN,width=10).grid(row=r,column=1)
    r = r + 1
    i = i + 1

root.title("GISSTUDYroot.title")
root.mainloop()

What is the problem?  Is there an alternative in Pythonaddins to display similar messageBox with an image?


texts:
I am following up on Handling data gap in new SoilGrids datasets which is a similar question about a previous version of SoilGrids250m v2.
There still appears to be gaps in northeastern Russia in the most recent SoilGrids250m version, including in the online viewer:
.
I have the same issue after downloading the entire dataset and re-projecting it to WGS84 or Winkel Tripel in ArcGIS.
Any clue how to solve it?


texts: Cannot export array bands in Google Earth Engine

texts:
I'm new to Google Earth Engine and Google AI Platform. A model was previously trained and published to AI platform and be imported to Earth Engine for image classification. The trained model successfully classified landcover in the specified area into 5 classes, however, I faced some problems when I was trying to export the classified images to my Google Drive.
Below is the script I have so far.
An error message of 'Cannot export array bands' jumped out when I tried to export the result.
How can I export the labeled image and the probability maps?
var model_1985 = ee.Model.fromAiPlatformPredictor({
  projectId : 'most-308413',
  modelName : 'Mostcnn_1985',
  version : 'v0',
  inputTileSize : [100, 100],
  proj : ee.Projection('EPSG:4326').atScale(30),
  fixInputProj : true,
  outputBands : {'landclass': {
  'type': ee.PixelType.float(),
  'dimensions': 1}
}
});

var predictions_1985 = model_1985.predictImage(array_image_1985);
var label_1985 = predictions_1985.arrayArgmax().arrayGet([0]).rename('label');
var probabilities_1985 = predictions_1985.arrayFlatten([["builtup","barren","cultivated","water","vegetation"]]);

Export.image.toDrive({
  image: predictions_1985,
  description: 'dl_classified_1985',
  folder: 'GEE',
  region: Clip_Extent,
  scale:20,
  crs: 'EPSG:4326',
})


texts: Running a script with QgsFeature from within the console

texts:
I am trying to run a script from the QGis Python console (similar to: http://spatialgalaxy.net/2012/01/27/qgis-running-scripts-in-the-python-console/)

The line feat=QgsFeature() operates fine from the Python console, however within the script it throws up the following error:

  Traceback (most recent call last):
    File "", line 1, in
    File "getLineFeaturesAndPoints.py", line 47, in write_polylines
      feat = QgsFeature()
  NameError: global name 'QgsFeature' is not defined

The code is saved in getLineFeaturesAndPoints.py file:

class WritePolylines:
def __init__(self, iface):
    """Initialize using the qgis.utils.iface
    object passed from the console.
    """
    self.iface = iface

def write_polylines(self, file_dir):

    # Get the active layer
    # qgis.utils.iface.mapCanvas()
    canvas = self.iface.mapCanvas()
    cLayer = canvas.currentLayer()
    #cLayer.name()

    # Get reference to the data provider
    provider = cLayer.dataProvider()

    # get all features (but no geometry)
    provider.select()

    # get attributes
    selectList = []
    columns = provider.fields()
    for key,value in columns.items():
        column = str(value.name())
        selectList.append(provider.fieldNameIndex(column))

    provider.select(selectList)
    features_file = open(path.join(file_dir, "features.txt"), "w")
    points_file = open(path.join(file_dir, "featurepoints.txt"), "w")
    writeFeature = csv.writer(features_file, quoting=csv.QUOTE_ALL)
    writePoint = csv.writer(points_file, quoting=csv.QUOTE_ALL)

    # get one feature id and geometry
    feat = QgsFeature()

Similar to the referenced HTML file, I invoke it as follows:

from getLineFeaturesAndPoints import WritePolylines
pl = WritePolylines(qgis.utils.iface)
myPath = 'C:\\Users\\Documents'
pl.write_polylines(myPath)


texts:
I have a shapefile which has the following PRJ

GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic latitude", NORTH], AXIS["Geodetic longitude", EAST], AUTHORITY["EPSG","4326"]]

As far as I understand, this means lat, long ordering. However, when I read the features using geotools, I get coordinates in long, lat ordering, at least that's what I see in my debugger after the features have been parsed. The problem is then that I do a reprojection to 4326 using the following code (I'm doing a generic shapefile import service):

        CoordinateReferenceSystem geo = CRS.decode("EPSG:4326", true);
        CoordinateReferenceSystem crs = feature.getFeatureType().getCoordinateReferenceSystem();
        MathTransform transform = null;
        if (!crs.equals(geo)) {
            try {
                transform = CRS.findMathTransform(crs, geo, true);
            } catch (FactoryException e) {
                throw new Exception("Unsupported projection type");
            }
        }

And this reverts all the lat long coordinates!

If it is correct to say that the axis order is important in the PRJ file, and assuming there is indeed an order mismatch between the PRJ and the SHP, how can I deal with this?


texts:
Using ArcMap 10.3, I want to digitized only the pale blue polygons in the raster into shapefile with "Raster to polygon" tool. This is the raster:

I get a shapefile with gridcode field, but i don't no how to isolate the right features of the pale blue color:

UPDATE:
this attribute table as 224000 features, i just show the first features in this picture:

I use unique values symbology and get this result (I deleted useless features from attribute table):

For clarity, I asked this question in https://geonet.esri.com/thread/165906.


texts: Export a configured Atlas with a Python script / command line?

texts:
I just want to export a configured Atlas (composer with Atlas option and cover layer enabled)

I tried this :

composerView = qgis.utils.iface.activeComposers()[0]
composition = composerView.composition()
atlas = composition.atlasComposition()
composition.setAtlasMode(QgsComposition.ExportAtlas)

atlas.beginRender()
for i in range(0,atlas.numFeatures()):
    atlas.prepareForFeature(i)
    composition.exportAsPDF('C:\xxxx\test.pdf' + str(format(i)))
atlas.endRender()

But nothing happens.

Have anybody a clue?


texts: Python processing script with clip in QGIS 3.2

texts:
I have been following the QGIS scripting templates, but have been unable to get vector processing to work with input and output layers. Working with QGIS 3.2.2.

The following script is my attempt at processing a clip. I am pretty sure the clip processes, but the 'Output layer' it produces is empty. Attempts at processing a buffer have also produced an empty vector layer.

from PyQt5.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
                   QgsFeatureSink,
                   QgsProcessingException,
                   QgsProcessingAlgorithm,
                   QgsProcessingParameterFeatureSource,
                   QgsProcessingParameterFeatureSink)
import processing

class ExampleProcessingAlgorithm(QgsProcessingAlgorithm):
    INPUT = 'INPUT'
    OUTPUT = 'OUTPUT'

    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return ExampleProcessingAlgorithm()

    def name(self):
        return 'clip'

    def displayName(self):
        return self.tr('Clip setup')

    def group(self):
        return self.tr('Example scripts')

    def groupId(self):
        return 'examplescripts'

    def shortHelpString(self):
        return self.tr("Clip test")

    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterFeatureSource(
            self.INPUT,
            self.tr('Input layer'),
            [QgsProcessing.TypeVectorPolygon]))

        self.addParameter(
        QgsProcessingParameterFeatureSink(
            self.OUTPUT,
            self.tr('Output layer')))

   def processAlgorithm(self, parameters, context, feedback):

       regveg = "Regulated Vegetation"
       source = self.parameterAsSource(parameters, self.INPUT, context)
       (sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,source.fields(), source.wkbType(), source.sourceCrs())

        cliplayer = processing.run("native:clip", {
        'INPUT' : regveg,
        'OUTPUT' : 'memory:clip',
        'OVERLAY' : parameters['INPUT']
    }, context=context, feedback=feedback)['OUTPUT']

        return {self.OUTPUT: cliplayer}


texts:
I have a topographic map of Washington as a TIFF-file from the US Geological Survey. And I'd like to clip the raster to the extent of the outer grid lines (78E/39N, 76E/39N, 76E/40N, 78E/38N). So that the map information (legend, projection, etc.) disappear. How can I realize that in ArcGIS?

Is there maybe a way to create a grid as a polygon and to cut the raster along this polygon? Which other possible ways are there?


texts: How to calculate the number of polygons along a line segment in QGIS

texts:
The image below is a representation of a digitised urban area derived from Google Satellite and OSM. The black lines between nodes are street segments and the irregular polygons are households.

I want to calculate the number of households along a street segment - the initial approach was to specify a buffer around each line segment in QGIS, and then use the 'Join attributes by location' to derive counts of households that fall within the buffer zone. However, I realise this is not a robust approach.

Is there an appropriate way of performing this calculation?


texts: Problem using polygonize / raster to vector in QGIS

texts:
I want to take my raster (multi ortophotos (.tiff) placed together using the "Build virtual Raster" option with the pyramids) and poligonize it. Something like this: gdal polygonize in python error, blank polygon created into an .shp file.

However, when I use the poligonize option an error appears saying that the process finished shortly after the beginning.

After having a closer look on my raster I found this:

Could this be the problem? How can I take them out of my layer?

After some research, I found the clip plugin but to use these I have to have a shapefile with the contour or something like that and I don’t have nothing. I must say I’m a new to QGIS and so to Geography :)

I thought about making my own contour, by drawing the boundaries, but could not found any plugin to do this. Is there a way to clean up my layer in order to pass it from raster to vector correctly? in shapefile so that I can open it, for example, in AutoCAD?


texts: I have four 15 arc second resolutions DEM files (15-C.tif, 15-D.tif, 15-I.tif, 15-J.tif) downloaded from viewfinderpanoramas.org. When I input them in QGIS I noticed that the 15-J.tif layer differs a bit from other layers in gray-scale intensity. Below is the image that illustrates the problem. I believe this lack of color matching can be overcome somehow. I would like to have adjacent DEM layers to match in color if possible.


texts: Joining shape to external data using ArcGIS API for JavaScript?

texts:
I'm developing an application that uses Shapes only to draw the map.
I would like to "Inject" the data at the time of viewing, since my data is dynamically changing.

Anyone know of any way of doing this, or if it is possible?


texts: Using st_intersects between tables in PostGIS?

texts:
I have two tables:

filter_polygon

and

myGeometries

In both tables I have geometry columns. myGeometries contains > 100 mln of rows.

In both tables I have geo indexes like:

CREATE INDEX filter_polygon_ix2
    ON filter_polygon.filter_polygon USING gist
    (geom)
    TABLESPACE pg_default;

Let's say that I have id of one polygon from filter_polygon table and would like to take all records from myGeometries table that intersect with it. My query looks like that so far:

select * from
filter_polygon.filter_polygon AS Filter,
myGeometries AS myGeoms
where
Filter.id = '0101'
AND
ST_Intersects(myGeoms.geometri, Filter.geom);

And the quey plan for it is:

'Nested Loop  (cost=261.66..25164.63 rows=1934865 width=251228)'
'  ->  Seq Scan on filter_polygon filter  (cost=0.00..5.24 rows=1 width=248951)'
'        Filter: ((origin_id_fieldvalue)::text = '0101'::text)'
'  ->  Bitmap Heap Scan on myGeometries  (cost=261.66..25139.51 rows=1989 width=2277)'
'        Recheck Cond: (geometri && filter.geom)'
'        Filter: (st_intersects(geometri, filter.geom))'
'        ->  Bitmap Index Scan on myGeometries_idx1  (cost=0.00..261.17 rows=5966 width=0)'
'              Index Cond: (geometri && filter.geom)'

The problem is that it takes like 20 min to get a result, sometimes even more when filter_polygon is complex and covers relatively big area.

Is there a way to improve that performance?

postgisVersion:

'POSTGIS="2.1.8 r13780" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.1, released 2013/08/26" LIBXML="2.9.1" LIBJSON="UNKNOWN" RASTER'

PostgreSQL Version:

'PostgreSQL 9.4.18 on x86_64-unknown-linux-gnu (Ubuntu 9.4.18-2.pgdg14.04+1), compiled by gcc (Ubuntu 4.8.4-2ubuntu1~14.04.4) 4.8.4, 64-bit'


texts:
I have a tiff image with 4 bands (Red, Green, Blue and Infrared) and I need to extract two of them.

I want to input one image and split it and make two images.
I want to do it with GeoTools since I use Java to develop.

How I do that?

PS: I can do it with gdal and command line. (e.g. gdal_translate -b 3 in.tiff out.tiff which 3 is the number of band)


texts: How can i combine the colours of one raster with values from another

texts: I have two Tiffs with data from the same area. One is in colour and when I use the identify shows a 'color index' number (eg. 141), then Color(a,r,g,b)(eg. 255,90,255,0). The other tiff contains a numerical value for the cell. I want a raster that contains both the colour info from the first tiff, and the numerical values from the second. Is there a way to add attributes from another raster or create a new raster that borrows from these two?


texts: How to convert from .SHP to .MAP for mapsforge?

texts:
We are using mapsforge, a Java component for Android.
This component use a proprietary format file (.map) that can be download from OpenStreetMap (OSM).

I need to use .shp maps, so I want to know how can I convert a .shp file into a .map file for mapsforge for Android


texts: Getting NDVI data only for March from 2000-2003 using Google Earth Engine?

texts:
I have a polygon and I want to calculate its mean NDVI values from 2000-2003 (and few more periods) but only for March
so the final result should be image which contains the mean NDVI values only for march between 2000-2003.

How can I do that?

This is the code I have written until now:

var dataset = ee.ImageCollection('MODIS/006/MOD13A1')
                  .filter(ee.Filter.date('2000-03-01', '2000-03-31'))
                  .select('NDVI')
                  .mean()
                  .clip(geometry)
                  ;
//var ndvi = dataset.select('NDVI')
var ndviParams = {min: -2000, max: 10000, palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ]};

 Map.centerObject(geometry, 7)
Map.addLayer(dataset,ndviParams,'NDVI2')


texts: Union multiple layers with OGRVRTUnionLayer

texts:
I'm trying to union overlay layers using OGRVRTUnionLayers. However, it keeps the new records on new rows as opposed to appending the data to the same row (So you would expect to get back all data but with polygons intersected from the different layers). I would like that if I click on one polygon, it returns all the attributes for that union overlay polygon. Here is my code which is placing the data on seperate rows for each layer:

<OGRVRTDataSource>
    <OGRVRTUnionLayer name="unionLayer">
        <OGRVRTLayer name="Layer1">
            <SrcDataSource>Layer1.shp</SrcDataSource>
            <SrcSQL>SELECT a FROM Layer1</SrcSQL>
        </OGRVRTLayer>
        <OGRVRTLayer name="Layer2">
            <SrcDataSource>Layer2.shp</SrcDataSource>
            <SrcSQL>SELECT b FROM Layer2</SrcSQL>
        </OGRVRTLayer>
    </OGRVRTUnionLayer>
</OGRVRTDataSource>

I would like it to do the same as the ArcGIS Union layer function.


texts: Line start and end nodes: determining what is a junction

texts:
I have a layer of points and polylines and want to determine which points are road junctions. I currently have all start and end nodes within 3m of my road network as a layer. However I cannot determine which ones are simple joins between sections or what is a junction, e.g:

In the image above, the highlighted point is not a junction, but it is where one section ends and another begins. Is there a way to distinguish between these and points like the others in the image without manual checking?


texts:
I want to download orthoimagery from this link:

https://geolibrary-maine.opendata.arcgis.com/datasets/595cf522a36b4c6593acc7a8d5db3190

But when i click download button, the result is;

How can i download high res ortho image?


texts:
How do i iterate through a shapefile using shapely. This is my code but it keeps outputting the same id. I think that it keeps looking at the same record.

fc = fiona.open("file1.shp")
print fc.schema
for counter in range(0,50):
    shapefile_record = fc.next()
    print shapefile_record['id']


texts: How to produce "SpatialPixelsDataFrame" from stacked RasterLayer?

texts:
How to produce "SpatialPixelsDataFrame" from staccked RasterLayer ?

I have a stacked rasterlayer and I want to convert them to SpatialPixelsDataFrame.

RasterLayer information is

> RAD.all[[1]]
class       : RasterLayer
dimensions  : 180, 240, 43200  (nrow, ncol, ncell)
resolution  : 0.01, 0.01  (x, y)
extent      : -123.305, -120.905, 36.705, 38.505  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : C:\Users\Michael\Documents\1998bil\1998ASC5min_ppt_spas1214_0001_19980202_0810_UTC.asc.bil
names       : X1998ASC5min_ppt_spas1214_0001_19980202_0810_UTC.asc


texts: Clicking on point to have additional points and geographical information show up in QGIS

texts:
I am mapping nonprofits in Costa Rica, using the headquarters as the main points on the map. I want end users to be able to click on the points and see additional information about each individual nonprofit like zone of impact (polygon), project locations (points) and additional offices (points).

To do this using QGIS, would I have to make a layer for each nonprofit?


texts: Good protocol for multiple conditional spatial and attribute based selections using ArcPy while avoiding saving to disc

texts:
I am running multiple conditional selections, both where clauses with attributes and location based selections.  I want to build functions, possibly within classes that are reusable and can pipe to each other or be scripted in different orders based on workflow.  I DON'T want to CopyFeatures_management constantly for intermediate features classes written to disc.  I'm hoping for Python type:

temp = load_data(fp_data)
temp2 = function1(temp)
temp3 = function2(temp2)
CopyFeatures_management(temp3)

Simple example of current workflow:

def select_by_location(fp_select, fp_location, distance_string, fp_out):

    arcpy.MakeFeatureLayer_management(fp_select, 'in_lyr')
    arcpy.MakeFeatureLayer_management(fp_location, 'location_lyr')
    arcpy.SelectLayerByLocation_management('in_lyr', 'within_a_distance', 'location_lyr', '5 Miles')
    arcpy.CopyFeatures_management('in_lyr', fp_out)

Call function with first feature layer:

select_by_location(fp1, fp2, '5 Miles')

Now take output and select with another feature layer:

fp_location2 = 'path/to/another/file'

arcpy.select_by_location(fp_out, fp_location2, '3 Miles', fp_out2)

What is a good protocol to avoid saving to disc constantly when in ArcPy world?


texts: 3D geological model from cross sections, preferably using Open Source software

texts:
I work with engineering geology, and my workflow includes making many cross sections of a given site, with data from maps, fieldwork (outcrops) and boreholes.

What I'm looking for is a way of taking many cross sections and making a 3D model, with open source software, even with lots of (CAD) drawing work, but preferably with few code or programming knowledge.

So, when I first noticed this software (GSI3D), that takes many cross sections (actually a fence diagram) and gives, as a result, a 3D model, made of solids and not just surfaces, I thought it would be the perfect software. The image below is from the software website.

Given that I don't have the resources to buy it, I was thinking if there wouldn't be a free alternative (open source software) to this approach.

So I think the approach would be:

Take the (CAD drawn) cross sections and combine them all in a GIS software, resulting in a fence diagram;
Interpolate between all of the cross sections, resulting in surfaces, or even better, solids, each one representing a geological layer or lithotype;
Make new cross sections between the former ones, to refine and adjust the model.

So, I think the 3rd step means that there is no need of a (very) complex interpolation algorithm, because the geological complexity would come from the cross sections, that are human made, not from the interpolation, that would be computer made.

The approach would be a first, simple, computer made interpolation, refined by manual adjustments.

To make it clearer, I'm not looking for a way of making a 3D geological model with complex interpolation algorithms to do all the work for me. I wouldn't mind having many iterations of cross sections to adjust and refine the 3D model, but I need a first (simple) interpolation so that I can work from it. The GSI3D approach seems to be it (I haven't used the software, but that's what I've understand).

I've read:

Looking for a free program for 3D visualization
Is there QGIS plugin to allow 3d visualisation of geological borehole data similar to functionality of Target for ArcGIS?
3D visualisation with open source GIS software
Combining Geological Grids and extruding between them?
a forum discussion
a paper (this one looks good, but the author didn't uses GIS)

among other GIS sources, and it appears that GRASS would be the best choice, but the discussion (when using GRASS) seems to be always about an algorithm to do the hard work or using an expensive, proprietary software.

If the answer is GRASS, what would be the tools to do it (as I sad earlier, is there a way of doing it with just a few programming knowledge)?


texts: Practical limitations of image classification of large rasters?

texts:
Been attempting (with no success so far) to implement the random forest classification in ArcMap on a reasonably large image. If I throw a raster with say 20-30 bands in it, works fine, just takes time. However my aim is really to get the tool to classify a 250-band raster, to really see how well it works on our data. Splitting the 250-band raster into multiple smaller band rasters defeats the purpose of what we're trying to achieve.

However this means our input raster (ESRI GRID stack format) is over 80Gb in size, and ArcGIS can't seem to handle it. So I am basically wondering if anyone has actually managed to implement any of the image classification tools on relatively large rasters? If not, is anyone aware of what the practical limitations of these tools are?

I have access to both Desktop 10.7 and Pro, run into problems with the large rasters on both and am running a desktop PC with decent specs.

Out of frustration, I have also tried to implement random forests using just a python script outside of the ESRI environment on a TIFF version of the GRID stack. Code works fine on small'ish inputs, but computer said "no" for the big files (the TIFF file is nearly 70 times the size of the GRID stack...).

Any advice? No access to a supercomputer at the moment, but if that's your advice, and you've successfully got it to work on one, happy to hear about it!


texts: QGIS 3.2: Processing toolbox tools from SAGA never run

texts:
I am having major issues in getting tools to run from the toolbox in version 3.2.  In essence all but very simple algorithms like 'Split RGB' return errors and do not run.  I cannot get any openCV or vigra tool to run. And I am not using the modeller.  Just opening the tool and setting params does not work.  The documentation seems to be in construction so it's hard to tell if I'm using this wrong or if there is a bug.

Before I start asking specific questions, do you need to install an OpenCV/vigra component in order to get these to run?  Do the procesing tools from external sources (Grass and Saga) actually work or are they experimental? A classic error message I get is:

Error: select a library.

Or

AttirbuteError: 'List' object has no attribute 'source'

EDIT:
I noticed that the QGIS 3.2 installer setup Saga 2.3.  But the OpenCV tools only appear in much later versions of SAGA.  If you open the SAGA version installed by QGIS, there are no openCV tools.  So I guess it's not surprising that when you run an OpenCV tool from the QGIS processing tools it cannot find its library. Is this worth an official bug report?


texts: Setting defaults in the 'Select By Attribute' SQL editor

texts:
Is there a way to set default values in the 'Select By Attribute' SQL editor (the where clause) ?

I want to set a permanent value that will appear every time I enter the 'Select by Attribute' window, like in this pic:

It doesn't matter to me if it'll appear in all layers that will be selected or only in specific ones.

Oh, and I'm using ArcGIS 9.3.1.

Thanks in advance.


texts:
How do I merge multiple binary rasters that each have different values and still retain the values in the output file?

Background Information:

I am working on improving unsupervised classification by applying 4 slope classes to it for a grad school project.For example unsupervised classification clustered playa together, but there are some playa pixels on the top of 4,000 foot high mountains. By applying the slope classes to it using raster calculator, I can isolate those pixels and change them to the granitic mountain class.

EDIT: I did it like this:
("unsupervised"==3) & ("slope">=30) - this output would become mountains ("unupervised"==3) & ("slope"<30) & ("slope">=5) - this would become alluvial slope. There are 9 unsupervised classes and 4 slope classes. 9x4=36 outputs.

This process produced 36 rasters using raster calculator. These all have a value of 0 and 1. I tried all the same things listed below to merge them, but all outputs just ended up with 0 and 1 in the value field except for the "combine tool". That makes sense to me.

So then I multiplied each one by a unique class number. Now I have 36 rasters that are reclassified, each one contains values of 0 and a corresponding class 1-8. So there are several files that are class 2. Several that are class 8 , and so on. They do not overlap that I know of. I want to get them all into one file with 8 rows, (or 9 rows with the 0 value).If I could even get 36 rows with their reclassed value retained and then use the "lookup tool" to dissolve them that would be ok too. I been googling for hours and have tried everything as follows.

Mosaic to new raster - results in 0 and 1.

Add to raster dataset and calculate statistics - did not work , don't recall why

combine - Can't do all at once, reassigns a new value. Yeah, I could make this work sort of, but no one who opens the file will know what is going on or how to symbolize it.

Raster Calculator - add them together "raster1" + "raster2" etc. Always have an error like something that was not in the input values for example, a row with a value of 16 when no input is over 8.

Raster Calculator - Merge([raster1],[raster2],[raster3]) Error says output raster empty or similar.

Append - The table looks better but the image is a disaster something is totally wrong.Seems to be too many 0 pixels, holes in the image all over.

I have VERY basic python skills and also have ERDAS, if those options offer a solution.

Here it is before I reclassed all the files that were only 0 and 1. This is how I want the final image to look, just in 1 file with 8 classes.


texts: Adding features to a vectorlayer using checkbox

texts:
I'm trying to add features to a vector layer using a checkbox.
    var featureOverlay = new.ol.FeatureOverlay;

$("input:checkbox").each(function(){
    var $this = $(this);

    if($this.is(":checked")){
        featureOverlay.addFeature(feature);
    }else{
        featureOverlay.removeFeature(feature);
    }
});

var polyCoords = [];
var coords = "95.61,38.60 95.22,37.98 95.60,37.66 94.97,37.65".split(' ');

for (var i in coords) {
  var c = coords[i].split(',');
  polyCoords.push(ol.proj.transform([parseFloat(c[0]), parseFloat(c[1])], 'EPSG:4326', 'EPSG:3857'));
}

var feature = new ol.Feature({
geometry: new ol.geom.Polygon([polyCoords])
})

var layer = new ol.layer.Vector({
source: new ol.source.Vector({
    features: [feature]
    })
});

var map = new ol.Map({
  target: 'map',
  layers: [
    new ol.layer.Tile({
      source: new ol.source.MapQuest({layer: 'sat'})
    }),
    layer
  ],
  view: new ol.View({
    center: ol.proj.transform([95.22, 37.98], 'EPSG:4326', 'EPSG:3857'),
    zoom: 4
  })
});


texts: Load vector features without adding the layer to the map in OpenLayers

texts:
Is it possible to make a vector source load and reload features without adding/removing its parent layer to the map? I need to do that after my viewparams url variable changes.

My vectore source looks as follows

var source = new ol.source.Vector({
    url: function(extent, resolution, projection) {
        var url  =
            'http://localhost:8080/geoserver/wfs?service=WFS' +
            '&version=2.0.0&request=GetFeature&typeNames=by:data' +
            '&outputFormat=GML3' +
            viewparams;
        return url;
    },
    format: new ol.format.WFS()
});

I've tried calling source.refresh() and source.changed() but it won't work.


texts: Mapnik/CartoCSS: make text labels overlap other layers (but not each other)

texts:
I'm using Mapnik 3.0.7 to render a basemap and a Shapefile layer that contains points. The labels for these points should overlap all the layers underneath, but not each other. Unfortunately, road labels from a lower layer are preventing some of them from being placed:

I've tried the following:

Adding allow-overlap="true" to the point TextSymbolizer, which makes the point labels overlap everything, including each other (bad).
Adding text-label-position-tolerance options to the underlying road labels, which I assumed would move them out of the way and allow the point labels to be placed. No effect.
Various text-comp-op options on the point layer and underlying layers, with no effect or undesired effects.

Is there any way to achieve this?


texts: How do you identify cuttings / crevices from a DTM or DSM using ArcGIS?

texts:
I am trying to identify where cuttings and crevices exist by using a DTM or DSM.

My ultimate aim is to find where roads and rail lines sit within a cutting i.e. the land either side of them is higher than the road or rail network.
If possible, I would like to identify the height difference between each road/rail and the top of the cutting, however, just finding the location of the cuttings would be great.

I have a line dataset and a polygon dataset showing me the roads and rails for my area.

My logic for this would be:

Find the crevices / cuttings from the terrain data and generate a raster with 0 and 1 values.  Where 1 is a cutting and 0 is everything else.
Vectorise this raster to generate polygons for the 1 values.
Union this with the road and rail dataset to find where the road and rail sit within the cutting.

By a cutting I mean something like this

Or this

I am not interested in this

It is at step 1 that I have hit a wall.

I would assume that some of the hydrology tools might be able to help me with this by modelling where water would flow. Does anybody have any suggestions?

I am using ArcGIS 10.2 with an Advanced Licence.


texts:
I'm using Google Earth Engine through the Python API and I would like to define my boundary area for the processing on the basis of a kml file I have. The kml basically describes the boundaries of a town, so it's an irregular polygon.
I was wondering if it was possible to run something like this:
MyVar=ee.ImageCollection("Sentinel").filterBounds(ee.Geometry.something("MyKml.kml"))

Otherwise I guess it is possible to "explode" the kml in a Polygon/MultiPolygon or something like that and pass that to ee.Geometry.
Any idea on how to do that?


texts: How to change programmatically the value of attribute in one of several feature in a vectorLayer

texts:
I have a vectorlayer that contains a lot of feature, I know that each feature has an unique Id

Each of them has some propertie like coordinateX, coordinateY, color, name

With a QTimer, each 2 sec I want to change the position of one of them  by change the coordinate using the QgsGeometry::fromPoint(x,y), i. e. the feature which has the Id=4 and update the Mapcanvas with the updated feature.

How Can i get the feature of a given Id ?
And
How can i proceed to update the mapcanvas with the change in the feature attributes ?

08/28 :added information
extract of my code , I am still working on it to find the solution :

//create point
QgsGeometry MyPoint;
MyPoint.fromPoint(QgsPoint(PosX,PosY));
QgsGeometryMap MyGeometryMap;
QgsFeatureId MyFeatureId=1;

MyGeometryMap.insert(MyFeatureId,MyPoint);

lDataProvider->changeGeometryValues(MyGeometryMap);
lLayer->updateExtents();

the data seems to be correct, PosX and PosY are random double, a qDebug show the value before store into the Geometry and it is correct
but the symbol don't move or if you prefer the symbol is not updated

Here is another code but still no change on the screen :

//get the feature of the specify Id
QgsFeatureIterator MyFeatIter;
MyFeatIter = lDataProvider->getFeatures();
QgsFeature MyFeature;
QgsFeatureId MyFeatureId;

while ( MyFeatIter.nextFeature(MyFeature)) {
    if (pTrackId == MyFeature.attribute("TrackId")) {
        qDebug() << "Feature Trouve  Id=" << MyFeature.id();
        MyFeatureId=MyFeature.id();
        MyFeature.setGeometry(MyPoint);
        MyFeature.setAttribute(Longitude,PosY );
        MyFeature.setAttribute(Latitude,PosX);
        MyFeature.setValid(true);
        qDebug()<< "TrackId=" << MyFeature.attribute("TrackId") << "  Label=" << MyFeature.attribute("Label")
                << "  PosX = " << MyFeature.attribute("PosX") << "  PosY=" << MyFeature.attribute("PosY");
    }
}

MyGeometryMap.insert(MyFeatureId,MyPoint);
bool lResult;
lResult=lDataProvider->changeGeometryValues(MyGeometryMap);

I try to set the geometry on the feature founded but no change,I've tried getFeatures on the dataprovider still the same result, coordinate change correctly but on the screen the feature don't change


texts: Can I drop my own templates into the OTB ExportWebMap Tool?

texts:
At ArcGIS for Server 10.1, there is a print task installed by default.

Looking at the web-help, it states:

  The PrintingTools service included with ArcGIS for Server references a
  preconfigured folder with 8 map layouts.  These include basic
  designs in the 8.5 by 11 inches, 11 by 17 inches, A3, and A4 paper
  sizes in both portrait and landscape orientations. If you want to use
  your own layouts, you can publish your own service for printing web
  maps

Despite this advice, I want to confirm if it is possible to include additional templates to this Out The Box service, as opposed to having to publish a separate service.

Looking at the Service in REST, I can see these templates.
I also found the corresponding templates under:
 C:\Program Files\ArcGIS\Server\Templates\ExportWebMapTemplates

If I drop a new 10.1 MXD into that folder, and restart the PrintingTools service (which now automatically clears the REST cache at 10.1), it will not pick up this new template.

Is there any way to configure the underlying GP service behind Export Web Map to accept new MXDs?

I did notice that if I edit one of these existing templates (e.g. add a North arrow) then these changes flow through.


texts: Why do 2 EPSG Codes exist for the same coordinate system?

texts:
I sometimes need to publish data in Gauss-Kruger zone 3 and was wondering why there are now the following EPSG Codes which define this CRS.

31467 - DHDN / 3-degree Gauss-Kruger zone 3
31463 - DHDN / 3-degree Gauss zone 3

The Parameters are identical. Then why 2 Codes?


texts: Using arcpy, update value in one field based on value in another

texts:
I am creating a level of service tool and I need to update the number of lanes field based on the operation type of the road. What I need it to do is: if the operation type equals 2W, then the number of lanes needs to be divided by 2.

for row in rows:
    for field in fields:
        if optype == "2W"
        arcpy.getValue(lanes / 2)


texts: How to get coordinate system info from shapefile using GDAL/OGR?

texts:
GDAL version 2.1.1

If I open a geospatial raster file, such as a GeoTIFF or DTED file, with GDAL, then I can get the coordinate system info for that file, as a string in OpenGIS WKT format.  For example:

GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]

The C++ code I use to get that info is simple and straightforward.

string path = "somewhere.dt1";
unsigned int openFlags = GA_ReadOnly;
GDALDataset * pDataset = static_cast<GDALDataset *>( GDALOpenEx(
    path.c_str(),
    openFlags,
    NULL, NULL, NULL ) );

string projWkt = pDataset->GetProjectionRef();

OGRSpatialReference srs( projWkt.c_str() );
srs.dumpReadable();

However, when I open a geospatial vector file, such as a shapefile, with GDAL, I get only a blank string for the coordinate system info.  The only difference in the C++ code I use is that I must also provide the GDAL_OF_VECTOR flag to GDALOpenEx().

string path = "elsewhere.shp";
unsigned int openFlags = GA_ReadOnly | GDAL_OF_VECTOR;
GDALDataset * pDataset = static_cast<GDALDataset *>( GDALOpenEx(
    path.c_str(),
    openFlags,
    NULL, NULL, NULL ) );

The shapefile I am opening does have a valid coordinate system.  I can confirm that by opening the shapefile in ArcMap or ArcCatalog.  Both of those apps are able to tell me the coordinate system and projection info.

So, why am I not able to fetch it via this GDAL API?

Is there a different mechanism for vector files?


texts: What kind of coordinate notation is this?

texts:
I've tried searching the web and a lot of questions here to find out which coordinate system this is:
[ 4.402223, 51.228429 ]
It should be near Antwerpen.
Does anyone recognise what coordinate system it is?


texts: Reading in different raster files with different bands using a loop

texts:
I've downloaded 6 multi-band raster files (each file has 4 bands) from the National Weather Service containing precipitation data for different time periods (30, 60, 90, 120, and 365 day periods). Here's where I downloaded the data from: https://water.weather.gov/precip/download.php

I want to write a loop to pull each of those files into R, with each raster band assigned to a different variable. Here's an example for just one raster file for the 30 day time period:

library(rgdal)
library(raster)

y <- format(Sys.Date()-1,"%Y")
m <- format(Sys.Date()-1,"%m")
d <- format(Sys.Date()-1,"%d")

NWS1_prev30 <- raster(paste0('/path/nws_precip_last30days_', y, m, d, '_conus.tif'), band = 1)
NWS2_prev30 <- raster(paste0('/path/nws_precip_last30days_', y, m, d, '_conus.tif'), band = 2)
NWS3_prev30 <- raster(paste0('/path/nws_precip_last30days_', y, m, d, '_conus.tif'), band = 3)
NWS4_prev30 <- raster(paste0('/path/nws_precip_last30days_', y, m, d, '_conus.tif'), band = 4)

I don't really want to copy and paste that 5 other times, so was hoping to write a loop. However, I'm receiving an error with my code as it is so far. Here's what I have:

pathname <- c(paste0("/path/nws_precip_last30days_", y, m, d, "_conus.tif"),paste0("/path/nws_precip_last60days_", y, m, d,"_conus.tif"),paste0("/path/nws_precip_last9days_", y, m, d,"_conus.tif"),paste0("path//nws_precip_last120days_", y, m, d,"_conus.tif"),
              paste0("/path/nws_precip_last180days_", y, m, d,"_conus.tif"), paste0("/path/nws_precip_last120days_", y, m, d,"_conus.tif"), paste0("/path/nws_precip_last365days_", y, m, d,"_conus.tif"))

# Bulk
for (j in 1:length(pathname)) {
  for (k in 1:4) {
    assign(paste0("NWS", k, j), raster(j, band = k))
  }
}

# Error in (function (classes, fdef, mtable)  :
#  unable to find an inherited method for function ‘raster’ for signature ‘"integer"’


texts: Creating small globe locator map in QGIS?

texts:
I've tried searching on gis.stackexchange and elsewhere but can't seem to figure this one out though it seems simple enough.

The type of locator I'm thinking of is similar to the one from the map in this post: Unattractive image artefact when exporting from QGIS's Print Composer

Does anyone know of a tutorial out there?


texts: Why output file of "mosaic to new raster" in ArcGIS is unreasonably large

texts:
I am trying to mosaic seven dem files together by using "mosaic to new raster". The largest one is about 100 GB, and the other six files are total 50 GB. I used "FIRST" mosaic method and the cell size is the same with the largest dem. And the boundary after mosaicing will not change, which is also the same with the largest one. So basiclly the output file should be about 150 GB or 200 GB something like that.

The mosaic work has taken a really long time, which is about 16 hrs and haven't finished yet. I checked the output folder and find there was only two files has been created, total 1.13 TB, which is unreasonable, because it is much larger than I expect. I am thinking if there is some temporary file there caused this and if the long time is reasonable for 100 GB DEM mosaic with 50 GB DEM?


texts: Select part of a polygon on one side of a line QGIS

texts:
I have a line shapefile of the outline of a water body. It was created by conversion of the polygon to line.

I am looking to add a layer of 500m standoff from the shore into the lake.

I have added a negative buffer of -500m, and that gives me the area which is accessible, however I want to just have that 500m as a polygon.

A picture speaks a 1000 words. See below.

What I am looking for is a polygon of the purple area.  When I buffer the line it buffers both sides, I just want one side.


texts: Buffer multiple files and distances on Python?

texts:
How does someone write a python script on arcpy.Buffer_analysis site package to buffer multiple shapefiles (e.g. roads.shp, schools.shp) with multiple distances (e.g. 500ft, 1000ft, 1500ft, 2000ft, 2500ft) and allowing them to have individual output file names?

New to python scripting.

Example:

arcpy.import

inputFeature = "roads.shp"

outputFeature = "roads_500ft"

bufferdist = '500 Feet'

arcpy.Buffer_analysis(inputFeature,outputFeature,bufferdist, "","", "ALL")


texts: Using the ArcGIS API for Python with a REST service to search content

texts:
I have a public REST service that has several layers and feature classes that I want to access. According to the documentation. Once I have the arcgis package imported, I can inject my outside REST service as such:
*Note: I am using the Jupyter Notebook to test this code.

from arcgis.gis import GIS
gis = GIS("https://hazards.fema.gov/gis/nfhl/rest/services/CSLF/Prelim_CSLF/MapServer")

Since this is a public service, I didn't have to include any username or password parameters. The service has several layers (item_type="Feature Layer) that I want to access; namely Special Flood Hazard Area Change layer, but perhaps others:

Using the example from the API docs, I tried to do a search in the service just to see all the layers and any info I might need in future calls. I am using the '.search' method from the content manager.

# search and list all feature layers in my contents
search_result = gis.content.search(query="", item_type="Feature Layer")
search_result

Every time I try this, I get an either a KeyError='num' or a TypeError: must be str, not int (if I try to add a query). What am I doing wrong here and what is the proper way to search within this REST service?


texts:
I have created a Reporting services report which has custom urls for  my CMV viewer based on the lat long for the feature in the record. I would like to click on it and open my CMV app a the custom Lat/Long and zoom level. Is that possible.

Here is an example:

http:///CMVApp/#/-122.71970/38.54753/19

The -122.71970/38.54753/19 would be what would change each time. I've tried this and sometimes it will zoom to the right location and other times it goes to the center of my map extent.


texts: QGIS timemanager NetCDF raster loads correctly but then run backwards

texts:
After timemanager loads correctly my time stamped NetCDF file (stamped in seconds using '%Y-%m-%d %H:%M:%S' format), the animation is played backward until it rewinds to the start.

I am using QGIS 3.4 and TimeManager 3.4.

This is not happening for every file. I could not figure out what is responsible of this behavior in my NetCDF file. Removing frame (outside QGIS) in a file that correctly loads can bring this faulty behavior.

Is there any known file format problem that cause this rewinding?

An example of a NetCDF file that does not load properly

netcdf knp14sku4_temp-04_2 {
dimensions:
    easting = 1104 ;
    northing = 1104 ;
    time = UNLIMITED ; // (62 currently)
variables:
    double easting(easting) ;
        easting:long_name = "easting UTM WGS84 UTM Zone -36" ;
        easting:standard_name = "easting" ;
        easting:units = "m" ;
    double northing(northing) ;
        northing:long_name = "northing UTM WGS84 UTM Zone -36" ;
        northing:standard_name = "northing" ;
        northing:units = "m" ;
    float temp(time, northing, easting) ;
        temp:_FillValue = -999.f ;
        temp:long_name = "lwir temperature" ;
        temp:standard_name = "lwirTemp" ;
        temp:units = "K" ;
        temp:grid_mapping = "transver_mercator" ;
    int time(time) ;
        time:long_name = "time" ;
        time:standard_name = "time" ;
        time:units = "seconds since 1970-1-1" ;
    char transver_mercator ;
        transver_mercator:spatial_ref = "PROJCS[\"UTM Zone 36, Southern Hemisphere\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",33],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",10000000],UNIT[\"Meter\",1]]" ;

// global attributes:
        :Conventions = "CF" ;
        :history = "Mon Jan 27 13:45:59 2020: ncks -d time,351,412 knp14sku4_temp.nc timeSplitFile/knp14sku4_temp-04_2.nc" ;
        :NCO = "4.7.2" ;
data:

 time = 1409060177, 1409060178, 1409060179, 1409060180, 1409060181,
    1409060182, 1409060183, 1409060184, 1409060185, 1409060186, 1409060187,
    1409060188, 1409060189, 1409060190, 1409060191, 1409060192, 1409060193,
    1409060194, 1409060195, 1409060196, 1409060197, 1409060198, 1409060199,
    1409060200, 1409060201, 1409060202, 1409060203, 1409060204, 1409060205,
    1409060206, 1409060207, 1409060208, 1409060209, 1409060210, 1409060212,
    1409060213, 1409060214, 1409060215, 1409060216, 1409060217, 1409060218,
    1409060219, 1409060220, 1409060221, 1409060222, 1409060223, 1409060224,
    1409060225, 1409060226, 1409060227, 1409060228, 1409060229, 1409060230,
    1409060231, 1409060232, 1409060233, 1409060234, 1409060235, 1409060236,
    1409060237, 1409060242, 1409060243 ;
}


texts: Google Earth Imagery not coinciding with OpenStreetMap in QGIS?

texts:
I am not sure what I am doing wrong, I did some research and Google Earth uses same Projection as OSM, but they differ. Whats even stranger, is that the layer added with quickmatchsearch from google earth does coincide with the OSM data, but the raster screenshot from the google earth app does not.

I exported the OSM data to DXF and need a raster image to put as background, but I need one that coincides.


texts:
I have a point feature class with residential information, but all that it has in terms of geographic information is the residence's house number.  There is no other information (street, zip, etc).  JUST house numbers.  The points fall more or less where they are located, but they are somewhat off.  I'm working off of a ArcMap basemap and I sometimes I can just move the points to areas that are obvious to me, for example if there's one point in a general area close to a house, and that's the only house, then its obvious to me that that point is that house.  However, i need a faster way to do this, and I need help with points that aren't so obvious, like for example when i have to assign a point to a residence and there are several residences to choose from.

I thought about geocoding, but I'm not sure how this would work when all I have to work with is a house number.

So I thought I could download a shapefile (maybe tiger data) with street names as well as zip codes, and do some kind of spatial join to where the points (house numbers) would be assigned the street name and zip code nearest the point.  and then create an address locator to match the points to the actual locations based on those attributes (street name, house number, and zip code).

Would this be a good way to tackle this? Does anyone have any other suggestions that might be better?  Could you maybe explain to me how I might do this? I'm pretty new to GIS and not very familiar with geocoding.
Also, can someone explain to me more or less how I should perform this spatial join?

Thanks!!


texts: OSM only label internal boundaries with Carto and Mapnik

texts:
Using the OSM workflow on switch2osm.org, we want to label all internal country boundaries but not external ones (coastlines).

There are 3 ways we could think of:

Is there any way to not show labels on polygon sides that have no neighbour (single sided)?
For example, instead of the usual OSM labelling:

Show something like this (done in ArcGIS with Maplex option to not allow single sided boundary labelling):

Using Carto, is there any way to position labels from one feature (e.g. boundary polygon) along another feature (e.g. internal boundary line)?
Use text-comp-op in Carto to only display the (non-existing) coastline label and ignore the lower level admin boundaries?

Or are there any other methods out there?


texts: QGIS Extract elevation data only for one country

texts:
I have a map built with a shapefile for country contours and a raster from .hgt files with elevation representation. How can I only select the area corresponding to the country?

I found some thing for matlab extract only the elevation of a country but not in QGIS. I'm thinking of using a polygon, but I am not familiar with how to do so. Or in a subtraction of the outside, but how to determine that? On the image, the red lines show the parts I want to eliminate.


texts: Viewing map tiles generated from Maperitive in leaflet or openlayers

texts:
for few days i have been working on developing my own osm tile server.
I tried it on linux as well as on windows,
no complete success till now.

I Have to make tile server on both windows and linux.

for windows i used Maperitive and succeeded in generating tiles and then i used
this guys code

http://gis.stackexchange.com/questions/84697/viewing-map-tiles-from-maperitive-in-leaflet-or-openlayers

now the prob is tiles are not showing up

Do i need to instal some Apache module for it ?
My tiles are kept here

C:\Profiles\Maperitive\Tiles

Where as i am opening html file from server

http://localhost:8080/mwork/maperitive_test1/index_leaf.html

which is kept in here

C:\Profiles\ApachePHP\apache\www\mwork\maperitive_test1


texts: pyQGIS changing Paletted Raster Layers to Grayscale

texts:
I've written a plugin which imports a number of colour geoTifs (map tiles) and works correctly but it would be great to display these as Grayscale after loading in.

The geoTifs load as paletted images (under Render Type) and I am looking to have them display as they do when you manually change the layer properties to set Grayscale to any of the three options (by lightness, by luminosity, by averaqe).

I've tried many things which haven't worked including things like layer.setDrawingStyle('PalettedSingleBandGray') which gives a very strange looking image - I assume it's only using information from a single colour channel?

I understand that these Grayscale options are relatively new to QGIS and hence their implementation through pyQGIS may not be straightforward, but I'm completely stuck with this now.

Following the response from xunilk, I've had a further play with this and whilst this does work, it doesn't quite do what I'm after.

Running xunilk's script for a colour tif gives:

Desired result (settings in image 3):

Setting grayscale OR desaturating the image as in image 3 give result like that in image 2.

Importantly, note that when SingleBandGrayRenderer is used (manually selected or through the script), even the manual selection of the grayscale option as in image 3 above gives a result that looks like image 1.

It seems that it needs to be loaded using QgsPalettedRasterRenderer with relevant rendering options specified?


texts:
I'm trying to automate a workflow that currently costs a co-worker about 45 minutes.  General run down is this person copies files from a network drive onto a local one, and process it through four separate ArcGIS models, the code below does not comprise all of the geoprocesses that are run against these files, but the point in which the bog down is happening.  As it stands, each model takes about 12 minutes to run.  The code below works, but it literally takes about an hour and a half on the same exact data.  I'm using arcpy.da.Walk because I figured I could save all the clicks and just run it in the background.  Pretty new to coding, but been doing it about 2 years - could it be the fnmatch?  My gut tells me it's the walk, and building the list - wondering if there is an alternative.

I have attempted to replicate the workflow (copying the directories locally, and not over the network).

def main():
    try:
        import arcpy, os, sys, fnmatch, traceback
        arcpy.env.overwriteOutput = True
        rootDir = arcpy. GetParameterAsText(0)
        outDir = arcpy.GetParameterAsText (1)
        mssn = '"12"'
        sensor = '"Geo"'
        dissolveFld = ["Date"]
        search = "*pan*.shp"
        outName = "test"

        shpMergeLst = []
        for root, dir, files in arcpy.da.Walk(rootDir):
            for filename in fnmatch.filter(files, search):
                shpMergeLst.append(os.path.join(root, filename))

        if shpMergeLst:
            outShp = os.path.join(outDir, os.path.basename(outName[0][:-4]) + "_Merge.shp")
            arcpy.Merge_management(shpMergeLst, outShp)

Don't have system specs, but currently running ArcGIS 10.1 Advanced


texts: How to get attribute of the leaflet marker?

texts:
is there a way to get and hold the option attribute of the leaflet marker?

marker = L.marker([lat,lng], {draggable:true,icon: L.AwesomeMarkers.icon({icon: '', prefix: 'fa', markerColor: 'darkblue',html:(100)}) });

In the above marker, how will I able to get the option attribute of ther marker like jquery .css method. such as, How to get if the marker draggable is true or false, or html value (100) like the .CSS attribute of jQuery $("#Div).css("color");


texts:
hi Im trying to access my raster file in QGIS by connecting postgresql database.
I can see my raster in qgis but I get this error :

Geçersiz Katman: GDAL provider Cannot open GDAL dataset PG: dbname=geoserver user=postgres port=5430 mode=2 schema=public column=rast table=boston: Couldn't establish a database connection Raster layer Provider is not valid (provider: gdal, URI: PG: dbname=geoserver user=postgres port=5430 mode=2 schema=public column=rast table=boston

as soon as I click it.

I'm using raster2pgsql for raster data and here is the sample of it:

raster2pgsql.exe -s 32637 -I -c -M -C -b 1,2 -M -l 2,4 C:\GDAL_Translate\boston_subset.tif -F -t 250x250 public.boston | psql -U postgres -d geoserver -h localhost -p 5430

Is there a problem in my raster2pgsql code ?
or
ıs there a problem occurs while connecting to QGIS ?


texts: Qt Designer and QGIS, using relation references

texts:
I have made my own form in Qt Designer and connected it to QGIS. One of my fields is a relation reference.
But the text doesn't fit in the combo box. You can see it in the picture.

I tried to change different settings in Qt Designer but it didn't help. Maybe someone can tell me how to change the size of a combo box in relation reference widget?


texts: Python ogr: why is GetFeatureCount() so slow?

texts:
I did some profiling for a program I wrote and discovered that ogr's layer.GetFeatureCount() is significantly slower than just using python's BIF len().

Here are my test results

def test():
        import ogr
        b_path = r'/../...'
        bounds =  r'%s/bounds' % b_path
        province =  glob.glob(r'%s/*.shp' % (bounds,))[0]
        driver = ogr.GetDriverByName('ESRI Shapefile')
        admin_shp = driver.Open(province,0)
        l = admin_shp.GetLayer(0)
        lc = l.GetFeatureCount() # first round
        lc = len(l)              # second round

if __name__ == '__main__':
    import timeit
    print(timeit.timeit("test()", setup="from __main__ import test"))

first round

228.255430937

second round

124.03532505

Can somebody explain why this is?


texts: Why is the QGIS print composer scale bar not right?

texts:
I went to re-do a map (of quite a large area of the globe) that I had already created, yet when I added the scale bar in the values were different and they definitely don't look correct, whereas before they were. The original map the scale bar went up to 8000km, which was correct, but now the same scale bar has the largest value of 1200km, when it should be 8000.

Unfortunately I overided the previous version of the exported PDF so I can't restore it.

I hope this makes sense.


texts: Understanding a typical output from a geotiff file

texts:
A typical output from geotiflib mey look like this:

Geotiff_Information:
    Version: 1
    Key_Revision: 1.0
    Tagged_Information:
       ModelTiepointTag (2,3):
          0                0                0
          1547250          5415950          0
       ModelPixelScaleTag (1,3):
          100              100              0
       End_Of_Tags.
    Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GeographicTypeGeoKey (Short,1): GCS_DHDN
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogEllipsoidGeoKey (Short,1): Ellipse_Bessel_1841
      GeogSemiMajorAxisGeoKey (Double,1): 6377397.16
      GeogSemiMinorAxisGeoKey (Double,1): 6356078.96
      GeogInvFlatteningGeoKey (Double,1): 299.152813
      ProjectedCSTypeGeoKey (Short,1): Unknown-31468
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      End_Of_Keys.
   End_Of_Geotiff.

PCS = 31468 (name unknown)
GCS: 4314/(unknown)
Ellipsoid: 7004/(unknown) (6377397.16,6356078.96)
Projection Linear Units: 9001/metre (1.000000m)

Corner Coordinates:
Upper Left    ( 1547250.000, 5415950.000)
Lower Left    ( 1547250.000,  941650.000)
Upper Right   ( 7186450.000, 5415950.000)
Lower Right   ( 7186450.000,  941650.000)
Center        ( 4366850.000, 3178800.000)

From what I understand ModelTiepointTag will tell me where in an actual map the pixels it identifies are located so in this case pixel 0,0 should be at 1547250, 5415950 (what values are these exactly?).

ModelPixelScaleTag will give me a relation about the images scale in real life. But what exactly do these values mean, and how do I extract a useful value from them (eg this distance in my image is these many meters in real life)?

Corner Coordinates is quite easy to understand, but again I dont know what do those values represent. Can I convert them to latitude longitude values?

If the above was too long, then basically what interests me is to be able to understand how can I correlate the values on the image with the actual world values they represent.


texts: Unsupervised classification in QGIS with Orfeo toolbox

texts:
I want to run an unsupervised classification of two landsat images in QGIS. Therefor I want to compute a k-means and a SOM (self organizing maps) classification. Unfortunately it doesn't work.

I tried to reinstall and install the new QGIS version (2.10) - and looked for Orfeo Toolbox Updates but unfortunately I couldn't run it (nor with different input data). Here is the logfile of the process in QGIS 2.10 (Windows 7)

 Algorithmus Unsupervised KMeans image classification startet...
C:\OSGeo4W\bin\otbcli_KMeansClassification
-in
"C:/Users/Manu/Documents/MASTERARBEIT/RemoteSensing_LandsatCDR/2010186_band4_BoaRef.tif"
-ram
128
-ts
100
-nc
5
-maxit
1000
-ct
0.0001
-out
"C:/Users/Manu/Documents/MASTERARBEIT/try/try.tif"
-outmeans
"C:/Users/Manu/Documents/MASTERARBEIT/try/try2"
Das System kann den angegebenen Pfad nicht finden.
Converting outputs
Lade Ergebnis Layer

The system can't find the path it states in the 3rd line from the bottom. (Path obviously is correct - and if defined error message stays the same.)
Same with the SOM algorithm in the following logfile extract:

    Algorithmus SOM Classification startet...
C:\OSGeo4W\bin\otbcli_SOMClassification
-in
"C:/Users/Manu/Documents/MASTERARBEIT/try/merge2.tif"
-tp
1
-ts
0
-sl
0
-sx
32
-sy
32
-nx
10
-ny
10
-ni
5
-bi
1
-bf
0.1
-iv
0
-ram
128
-rand
0
-out
"C:\Users\Manu\AppData\Local\Temp\processing\9d88158a0c1d42aeaf0ffa70e4744d0d\out.tif"
-som
"C:\Users\Manu\AppData\Local\Temp\processing\8bb385e883cb44039abe7125bb2bbd35\som.tif"
Das System kann den angegebenen Pfad nicht finden.
Converting outputs
Lade Ergebnis Layer

I suppose the toolbox doesn't work properly as I tried different input data, different QGIS versions, different output paths (from above are temporary; but I also tried to save the output directly on disc - ergo providing a path),...

Should I contact OrfeoToolbox directly?


texts: Using Python to publish a JavaScript webpage

texts:
I am building a series of HTML web maps for a client (5000 location maps all HTML). I have a python script to execute a basic set of HTML tags to build a specific web page (Call this page 1). The script is simple, upon running it will return my web page. Now, I am trying to do a bulk string find/replace to modify the HTML elements (Different characteristics for each map) I want python to search for a specific string in my code, then look up data from a spreadsheet to insert the HTML elements into it. But not sure how to go about this, I have a working python script that populates an HTML file. Next step is the look up and find/replace, do I use JavaScript or python for this and how do I go about it?

import os, sys
import webbrowser

f = open('XXX.html', 'w')
message = """<!DOCTYPE html>

    <title> Web Map Viewer</title> <button onclick="printMe()" style="float:right">Print Maps</button>
    <meta charset="utf-8" />

Location Web Viewer Map

    <fieldset >
        <legend>XXX Location Snapshot:</legend>
            Status: Existing | Latitude: 42.557, Longitude: -94.680, Elevation: 7348 ft. | Address: 500 State Highway 156 | Contact Phone: 123-345-3333
    </fieldset>

Map 1 - Substation and Nearest Vicinity:

Map 2 - Regional Surrounding Area / Closest Major City:

"""
f.write(message)
f.close()
webbrowser.open_new_tab('XXX.html')


texts:
1) I tried to define a custom CRS in QGIS without success. I fill all fields and try to change project CRS. (all parameters is good) But I can't find custom CRS in list. I think that this is due to this error.

2) This window opens at startup of the program. I tried to install all versions of QGIS from 2.6 to 2.18 and this error didn't disappear.
What am I doing wrong?

3) I remove qgis.db and the error message and the error message don't open now. But I still can't define new custom CRS.I fill all fields with this parameters.

When I try to change the project CRC I can't find it in list.


texts: How to measure and mark a line or point as a template QGIS?

texts:
I'm digitizing some of our sewer assets and need to be able to measure distance and create points.

So, for instance, I have a sewer main that is 65m long. I need to mark a point at the 32m mark, then a further 6.2m.

The main is represented by a polyline and already exists. I want to be able to create the points in a separate layer. In my innocence, I'm imagining something like a ruler tool where I can click to drop a point at the distance locations I need.

I have read the examples where a regular pattern of points are created, and have tried using the v.to.points tool in GRASS, but none of these give me the result that I want.

I can do something like this in MapInfo, but would like to see if I can do it in QGIS.


texts:
I've successfully loaded a bunch of shapefiles into a PostGIS database. I would like to be able to retrieve specific geometry and layers as dxf files based on SQL queries.

I'm using PostGIS, GDAL, PostgreSQL, and Python from the OSGeo4W suite on Windows.

I'm imagining that I should be using ogr2ogr a bit like this:

$ ogr2ogr -sql "SELECT * FROM some_table" -f "dxf" outputfile.dxf -overwrite "PG:dbname=mydb user=dbusername password=dbpassword"

but when I do that with the appropriate values, I get a series of errors like this:

ERROR 1: DXF layer does not support arbitrary field creation, field 'objectid' not created.

Resulting in an error for every attribute column in the layer some_table But I do get a dxf file.

If I've already created a .dxf file, and I rerun the command above (including the -overwrite option), I get this error instead:

Layer sql_statement not found, and CreateLayer not supported by driver. ERROR 1: Terminating translation prematurely after failed translation from sql statement.

I have already checked to make sure that the GDAL_DATA environmental variable is set correctly. (to C:\OSGeo4W\share\gdal)

It's certainly possible that I'm approaching this task incorrectly because I'm new to all these tools, so I'm open to any suggestions and tips about how to make .dxf files from GIS data stored in PostgreSQL based on SQL queries.