texts: Using color ramps (cpt examples and wikicarto 2.0) in QGIS 2.4?

texts:
I´ve been trying to use some of the color ramps you can find at cpt-city but I have problems with the value ranges (right word?) for each color when loading them. I just don´t understand how to get the settings in the Style window right so the desired colors appear as the authors intended.

I´ll give an example to show my problem:

I created the colorscheme “Wikicarto 2.0” as a qml-file by putting the example code into a texteditor and saving as .qml as suggested here:
http://en.wikipedia.org/wiki/Wikipedia:Graphics_Lab/Resources/QGIS/Create_a_topographic_background#Existing_Wikipedia_color_palettes

Now if I load this qml as my style in the properties window of my DEM everything looks fine, especially the Values used for each color. :)

If I change settings in the part where it says “Generate new color ramp” or “Load Min/Max values” it only throws everything out of order, the values for the colors change and I don´t know what I would have to change in these settings to get them back to the “correct” values. I can only load the qml again to solve the problem.

Now I also want to use other colorschemes like the “DEM print” which can be found here:
http://soliton.vm.bytemark.co.uk/pub/cpt-city/td/tn/DEM_print.png.index.html

I guess the values for each color in this scheme are specified too, just like in “wikicarto 2”…?

I downloaded this scheme as svg file and imported it into QGIS 2.5 with the plugin “Svg2color”, then exported it as a xml  “style” file.
This xml I then imported in the Style Manager.

I can now select the colorramp from inside the properties window of my DEM under “Generate new color map” in the dropdown list.

If I do this I then have to “Load values” in the section below and also click “classify” to get the style to show by clicking “Apply”. So somehow there is no “set” Information for value ranges inside this imported colorramp from cpt.

I also tried to just rename the xml file to qml but then I can´t load them through “Load Style”in the properties window (Error).

I then tried to use the “Wikicarto 2.0” as it can be found on cptcity in order to better understand how this all works by using the colorscheme I know from a different source and fileformat. On cpt city the “wikicarto 2.0” can be found here:
http://soliton.vm.bytemark.co.uk/pub/cpt-city/wkp/template/tn/wiki-2.0.png.index.html

If I import it into QGIS via the downloaded svg file and save it as an xml I again get the problem with the “correct” value ranges not being specified when selecting the WikiCarto 2.0 colorramp from the dropdown menu under “Generate new color map”. These values depend on how I set the Min/Max Values, classes and so on and only appear when I click on “classify”.

Maybe I somehow have to incorporate the info like segment count and valueranges given below the preview of the color ramp on cpt-city like this one here? I just don´t know where to enter it :(

To summarize:
I thought that inside the prepared colorramps I found there is also the information about the value ranges (corresponding to heights) for each color as desired by the author and I just want to know how to get this information from the downloaded svg or xml file into a qml-file so I can load it as easily as the wikicarto 2.0 scheme made by putting the code from the wikipage into a qml file.


texts: Extract values from one raster based on a conditional function in another raster

texts:
I am relatively new to Arc and having a bit of trouble with Extraction. I am using ArcMap 10.4.1 and I have two raster files: one has population data (density per 1km cell) and the other contains a specific disease risk (scale 0-1). Both rasters are the same size (size of cells, number of columns and rows). I want to be able to extract the values from the population raster based on certain values from the disease risk raster; ie how many people reside in areas of 0.6 to 0.7 risk of disease?

I have had a look at the different Extraction options but none offer this option (I think). Do I need to first extract disease risk by attribute using a SQL expression (e.g 0.6-0.7)->transform to polygon->extract population values by 0.6-0.7 risk polygon? Seem like a lot of steps so wondering if there is a more direct way?


texts: Creating Standard Deviation Index (SDI) cell by cell in DEM

texts: I am fairly new to QGIS and have been struggling to find a solution. To put it simply I am trying to find a way to display the Standard Deviation Index on individual cell basis by comparing it to its surrounding 8 cells. The purpose is to help narrow down possible habitat of a bird species. How can I achieve it?


texts:
I have 2 layers on PostGIS: a point layer named 'tdia' and a polygon layer named 'com':
-the com layer contains the city name. (Polygon)
-the tdia layer contains bus stops. (point)
I want to create a trigger on my layer tdia that allows to recover the name of the city after creating a bus stops.

My code does not work:

CREATE TABLE tdia (gid SERIAL NOT NULL PRIMARY KEY, numero
VARCHAR(5), CODE_INSEE VARCHAR(5), NOM_COM VARCHAR(50), date DATE, geom geometry(point, 2154))

CREATE OR REPLACE FUNCTION maj_tdia()
RETURNS "trigger" AS
$BODY$
BEGIN
UPDATE tdia SET CODE_INSEE =
(SELECT com.CODE_INSEE
FROM com, tdia
WHERE st_intersects (tdia.geometry, com.geometry));
END;
$BODY$
LANGUAGE plpgsql;

CREATE TRIGGER maj_tdia BEFORE INSERT OR UPDATE
ON tdia
FOR EACH ROW
EXECUTE PROCEDURE maj_tdia();

What am I doing wrong?


texts: How to deal with versioning in OpenStreetMap?

texts:
The topic of managing geospatial data in a more general sense has come up before here. The topic of versioning was mentioned there as well, but not really dealt with.

Traditional geospatial data collection and maintenance only need to deal with versioning internally, as the database is only updated from within the organization. This is not the case in crowdsourced geodatabases like OpenStreetMap. There, anyone can come along and add, modify or delete objects. In OpenStreetMap this is dealt with in a rudimentary way: each object has an integer version number, and only the object with the highest version is exposed in the live database. The database uses optimistic locking, so users must resolve all conflicts that occur when uploading contributions manually.

This all works reasonably well as long as human contributions through the editors (JOSM, Potlatch) are the only mode of contribution - but they aren't. Increasingly, imports of open public sector data are conducted. These make for more complex versioning issues. Consider the following scenario:

A building object is being imported from an open public sector dataset
The building receives some modifications by human contributors (attributes, geometry, or both)
A new version of the public sector data becomes available and is imported.

Currently, in step 3. the human contributions would be lost, unless each building that received community modifications is manually merged with the new import.

How can OpenStreetMap deal with this situation? Do we need to look at distributed version control in software development? How can methods of DVC be adapted to deal with distributed spatial data maintenance?


texts:
How do i get a new line in a Geoserver style CSS ?

{
label: [name] (need new line here) [pop];
}

Details: [name],[pop] are attributes of the layer

Is there a comprehensive CSS documentation with examples for Geoserver? I want the label to look like this


texts: Why does clipping polygons result in no features using QGIS?

texts:
I'm using QGIS 2.6.

I'm trying to clip the polygons in one layer by another but cannot get the results I expect. I have a layer with wetland areas, and another layer with the boundary of my client's property. I want a layer showing only those wetland areas that lie within my client's property. If I run the 'Clip' tool, the resulting layer has nothing on it. It makes no difference which layer I use as the input layer and which as the clip layer, the result is always no features.

Could anyone explain what's going on and what tool might do this operation if not 'Clip'?


texts: Is it possible to alter field display precision whilst using a coded domain in ArcGIS 10.2?

texts:
Firstly a bit of background:
I'm currently using fields of type double to hold the results of our chemical analysis. Different samples were analysed for different analytes and several values were reported below the limit of detection or not analysed. In an attempt to make this clear I have used a coded value domain to swap specific values to text as in the image below. I know this isn't strictly what the coded domain was intended for.

The problem:
What I want to be able to do is change the attributes field display to a specific number of decimal places. i.e. I want the numbers in the TotalNitrogen column to display to 3 decimal places.
What I have tried so far:
I have used numeric option under the field properties to specify the number of decimal places but this doesn't appear to work.
Additional caveats:
As I do not know what the end user of the dataset is likely to need I don't want this precision change permanent. I'm not adverse to using a different solution than using a mixed coded domain environment to qualify results (i.e. specify value is Below Limit of Detection) if a solution exists.
I am using ArcGIS Basic Licence level version 10.2.2.3552 in a Windows 7 environment if that helps.
Additional information requested in comments.
@Dan Jurgella The image below is how the coded domain is set up.


texts: Service Area (ArcGIS), including slope as a cost

texts:
I'm trying to make a service area network analysis containing the time that is needed using a bike starting from one point.
The attached map shows this service area when the bicycle speed is 20km/t, and a travel time interval of 5, 10, 15, 20 and 30 minutes.

But I want to include slope as a cost in the analysis, and calculate a speed based on the slope.
The speed will then be used to calculate a travel time for each segment.

The slope/speed calculation can be something like this (interpolation between):

downhill steeper than minus 8% = 40 km/t

-4% = 30 km/t

0 % = 20 km/t

4% = 10 km/t

upphill steeper than 8% = no bicycle

My biggest issue is that I don't know how to define the road links' direction.

What I've done so far:

1) Import road network

2) Dissolve road network, and split network into sections of 50m (since I will calculate the average slope by intervall of 50 metres)

3) Extract start and end point of each section (will be used to define altitude data for startpoint and endpoint)

4) Get xy-coordinate for each point

5) Get elevation data for each point from DEM model (metres above sea level)

6) Copy elevation data and xy-coordinate from each point to the 50m lines

I now have a network consisting of road segments of 50 metres with xy-coordinates and altitude data for each start and end point.
I now want to calculate the slope based on (elevation end-point minus elevation start-point) divided by segment lenght.
This will give me the slope (average), BUT I do not know the direction (upphill/downhill).

I thought about using the coordinates, but the road can curve a lot so based on xy-coordinates I actually do not know if I am travelling from or towards the starting point..

Any idea how to define the direction of the roads away starting from one point?

For example, is there a method to sequentially number nodes/points starting from one point?


texts: Is there any 'Auto-complete polygon' tool in QGIS?

texts: Is there any 'Auto-complete polygon' tool in QGIS, similar to the one existing in ArcGIS?


texts: Correctly filling in Parameters of QGIS Toolbox Tool "selectbyexpression"?

texts:
I'm trying to use the tool "qgis:selectbyexpression" in one of my own scripts.  But it doesn't work.  My question is simple: "Why?"  :)

My hunch is with the 2nd parameter, the 'expression' (see "HOW I GOT HERE" below).

MY CODE:

Input = QgsMapLayerRegistry.instance().mapLayersByName('Streets')[0]
iface.setActiveLayer(vl)

processing.runalg('qgis:selectbyexpression', Input, '"STRNAME" = streetslist[k]', 1)

Input = 'Streets', a shapefile with columns of civic numbers, street names, etc.

streetslist = a list I made with a few of the street names that I know match some of the streets in the "STRNAME" column of 'Streets'

HOW I GOT HERE

As with every Q tool from the Toolbox that's new to me I open it up and read up on it in the "Help" tab.  Personally, I often find these Help tabs rather unhelpful.  Case in point:

It shows the correct order of parameters but without explicit examples of how to fill in them in I'm left guessing.

So I reasoned that the "selectbyexpression" Toolbox tool would work the same as the GUI tool of the same name (top of the Q Project Window, in the Attributes Toolbar?  I'm sure you're familiar).

When you run the GUI tool expressions are built something like this:

I assumed the entire 'expression' needs to be surrounded by single quotes because I've seen other tools run that way.  Double quotes around the column name and single quotes around the value, right?  So I tried the following:

1) '"STRNAME" = 'Main''

This throws the following error:

File "<input>", line 1
    processing.runalg('qgis:selectbyexpression', Input, ' "STRNAME" = 'Hollis' ', 1)
                                               ^

2) ' "STRNAME" = 'Main' '

Exact same error message as previous except oddly the ^ appears to be moved to the right so that rather than pointing at the comma before "Input" it's pointing at the white space before "Input".

I even tried double quotes around the value (the street name) instead of the single quotes:

3) '"STRNAME" = "Main"'

Interestingly, this one runs without error but nothing actually gets selected.

I could go on trying various permutations of single quotes, double quotes, spaces, etc. but I figured the answer is quite simple and someone out there will be kind enough to share.


texts: How do I get around license error when creating custom tool?

texts:
I created a custom geoprocessing function tool according to this link . . http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#/Building_a_custom_geoprocessing_function_tool/00010000049w000000/

I went through the following instructions in order to create a custom tool in a toolbox, but got the error "You do not have the necessary license to execute the selected tool" when attempting to run the custom tool . . .

1) Register the DLL (in the Debug folder) with ESRIRegAsm utility. ESRIRegAsm.exe is located at, ~:\Program Files\Common Files\ArcGIS\bin folder. Use this example . . . Do the following to register the assembly:
               Desktop—ESRIRegAsm /p:desktop "C:\Test\MyAssembly.dll"

           Engine—ESRIRegAsm /p:engine "C:\Test\MyAssembly.dll"

2) Create a custom toolbox in a folder or a geodatabase.
ArcGIS for Desktop users have the option to add the tool to a custom toolbox manually in ArcCatalog, ArcMap, or ArcGlobe. Right-click a custom toolbox, choose Add, click Tool, expand my function factory, select my tool, and click OK. The tool is added to the custom toolbox.


texts: inputing road network as nodes and arcs in java

texts:
I want to input a road network as nodes and arcs. Some of the nodes are intersections, others serve to describe the shape of the road. Road network should be modeled as an embedded graph consisting of nodes and arcs.

I want to model the input road network as a graph G with node set N and arc set A, and with the following information:

Each node in the set N is a point with geographic coordinates. Each node has references to the incident arcs of the graph.
Each arc in A has some importance, representing the type of road (dirt road, ..., highway), and it has references to the two incident nodes.

Please help me with the code in java for doing the same.


texts: PostGIS importing shapefile error w/ official manager

texts:
I'm having some trouble when trying to import an shapefile of the Sao Paulo Brazil's state (I downloaded it from the web and I'm using the own PostGis manager to import it).

Here's an image!

The connection is already on.

In the "Options" menu, i had to change from "UTF8" to LATIN1, cause if I didn't, I'd get this error message: "Unable to convert field name to UTF-8 (iconv reports "Illegal byte sequence"). Current encoding is "UTF-8". Try "LATIN1" (Western European), or one of the values described at http://www.gnu.org/software/libiconv/."

I also wanted to now what is SRID and wich value should i put there.

When I changed to LATIN1, I get this other error message: "Failed in pgui_exec(): ERRO:  sequência de bytes é inválida para codificação "UTF8": 0xe3 0xa7 0xe3". In english, it says the byte sequence is not valid to UTF8 codification.

I'm a bit desperate cause I've never used GIS before.


texts: I have a raster layer (coloured squares) and a vector layer that represents water bodies (blue coloured polygones; see image below). I am looking for a way to exclude / remove the cells that are covered by more than a certain percentage of water (say 60%). How can I do it using gdal and / or python?


texts:
I have problems to release my gdb in ArcGIS10. My tool works with an local gdb, but I failed to release it after the process.
In Ag9.3 I can release it simply by Marshal.ReleaseComObject(workspace), but this doesn't work in AG10. Also releasing all FeatureCursors doesn't help. Only closing the whole application "deletes" the lock.
Does anyone know how to release my gdb and "delete" the lock without closing the application?


texts: CARTO.js (V4) - How to style different values of a polygon

texts:
Using Carto.js V4, I have added a layer that shows fire risk levels (low, medium, high and extreme) for South Africa.

This is the code I have used to call the layer from my Carto library:

  const veldfireDataset = new carto.source.Dataset(
          veldfire_risk_2010
  );
        const veldfireStyle = new carto.style.CartoCSS(
          #layer {
          polygon-fill: #162945;
            polygon-opacity: 0.5;
            ::outline {
              line-width: 1;
              line-color: #FFFFFF;
              line-opacity: 0.5;
            }
          }
  );
        const veldfire = new carto.layer.Layer(veldfireDataset, veldfireStyle);

How do I fill the polygon depending on the "overall" column's value? And ultimately also add a legend to explain what each colour relates to.


texts:
I have millions of points having temperature and humidity information covering the Middle East and North Africa (MENA).

How can I create animated heatmaps/hotspot maps of change over time? The data have a time and date field. I would like to do this using QGIS, either version 3.2.2 or 2.18.

An example of what I am trying to accomplish: https://www.climate.gov/maps-data.


texts:
I have a Multipolygon-Layer stored in a PostgreSQL/PostGIS database.
The primary key field 'gid' is SERIAL.
When trying to add a new polygon in QGIS (LYON 2.12.1 on OSX El Capitan), the attributes form shows me:

nextval('test_gid_seq'::regclass) for the gid field (which is the correct sequence to autoincrement the field).

When I try to save the edits, I always get an error:

ERROR:  duplicate key value violates unique constraint "test_pkey"
    DETAIL:  Key (gid)=(1) already exists.

This is strange, since the actual value in the sequence is 114. Why does QGIS tries to use 1?
Is this a bug in QGIS? Or does anybody know how to avoid these problems?


texts: Setting default value for Python script tool properties using ArcPy?

texts:
I have a model folder which contains several shapefile, but only one named oil.shp. Is there a way to automatically set the default value for the Add Oil Layer based on the model folder and the shapefile's name? It can be done by using os.path.join(model_folder, 'oil.shp') if I don't set Add Oil Layer as one parameter. However, I want to leave the parameter there in case someone else wants to play around with another oil.shp as needed. The Script Properties comes from the Script tool.

I am using ArcGIS Desktop 10.3.1 Advanced License.


texts: How can I extract the coordinates from one layer and import them in a diferent layer when key field is the same?

texts:
I'm working with 2 different layers.
   First one is a csv file with 3 fields. Name, longitude and latitude.
   Second is a text file with the field name.

I want to assign the longitude and latitude fields to the second layer when the name coincides, anyone knows how to do it?


texts: Finding UK population GIS data by city for free?

texts:
I am looking for some gis data for free of the UK population by city.

I have found some in excel files but I am looking for other data formats(shapefile? CSV? I am not sure), data divided in smaller administrative districts.


texts: Can you use active layer as input for dynamic text?

texts:
I'm wondering if you can lift the name of a layer in the ToC (within a group) and use it in dynamic text.

Is there existing syntax to do this?
Perhaps ...

<dyn type="ToC" property="activeLayer"/>

Using ArcMap 10 sp3.


texts: Where can I find documentation on the scaling levels of Open Railway Map tiles

texts: I am using the Open Railway Map xyz tiles (https://a.tiles.openrailwaymap.org/standard/{z}/{x}/{y}.png). The level of detail changes at various zoom levels. Where can I learn what the defined levels of detail are at each zoom level?


texts: How to handle service roads off of highways?

texts:
I am setting up a road network for a our rural municipality in Canada, and I am wondering how to handle the service roads that run parallel to our highways. The buildings along these service roads have addressing based on the highway, so it wouldn't make sense to create names for these roads.

I am wondering what other people do in situations like this, and if there is a "best practices" way to handle this. None of them have official names, so would naming the roads (simply for presentation purposes) something like "HWY 52 Service Road" be a good idea, or simply leave them unnamed, but still have them drawn out so that they appear in maps.


texts:
I want to split a QgsPolyline geometry into two separate lines at a given QgsPoint (by projection onto the line), adding the point to the line as end point of the incoming line and start point of the outgoing one.

I got the following code that works for 2D geometries (creating a QgsPointXY point and a QgsPolylineXY poly from the input QgsFeature):

def breaklineatpoint(pointfeat, polyfeat):
    point = pointfeat.geometry().asPoint()
    poly = polyfeat.geometry().asPolyline()
    geo = QgsGeometry().fromPolylineXY(poly)

    (sqdist, projpoint, aftervertex, leftOf) = geo.closestSegmentWithContext(point)
    newline1 = poly[:aftervertex]  # from start to break point
    newline1.append(projpoint)  # append point
    newline2 = poly[aftervertex:]  # from break point to end
    newline2.insert(0, projpoint)  # insert at position 0

    return newline1, newline2

However, it drops the feature's Z and M values. How do I translate this function so that it can handle 3D geometries (QgsPolyline instead of QgsPolylineXY and QgsPoint instead of QgsPointXY)?

I already fail creating a QgsPolyline from a QgsGeometry. Using a QgsLineString by QgsGeometry().constGet() is no alternative as I then can't apply the Python list cutting and slicing functions.


texts: QGIS Calculate distance between polygon vertices / find vertices that are within distance

texts:
I have polygon layer and i need to find (select, so i can export selected to new layer) all the vertices that are in provided distance to each other (like 2 meters).

I already tried:
1. Export polygon to point (Extract Node)
2. Run delete duplicates (as the extract node creates duplicates).
3. Create distance matrix (this one is problematic as it shows distance in datasource units which is dd and i want it to be meters)
4. Join tables and then filter them using distance < value

Is there a better way to achieve this, like a tool?


texts: How to have QGIS show the number of points that overlay the raster polygon?

texts: I have a number of single-value raster layers defining a area mask and a point layer pulled from PostGIS. Is there a way to have Quantum GIS show me the number of points that overlay the raster polygon? Or is my only option to convert the raster to a vector, pull into the database, and run a query against it?


texts: How to open a shapefile using Python in QGIS?

texts:
I have tried to open my shape files with the help of a tutorial (http://spatialgalaxy.net/2012/01/27/qgis-running-scripts-in-the-python-console/) but unfortunately I couldn't open my shapefiles for my first work in QGIS. I put the directory of the shapes and the script in my python path how the tutorial recommended. So above are two pictures showing the problem by script and in command prompt:

Are there any other solutions (in a form of python language) to open shapefiles? I am having headeches because when I'm trying to do like the website exposes I couldn't open the shaoes and doing the commands of importing Loader I got messages of errors. Hope someone here help me!

///////////////////////////////////////////////////////////////////////////////////////////

I have tried this afertnoon reproduce all the process dicussed here. I put the the script in the folder of my python path I think is correct. As I said any times later when I use the command "sys.path" many directories appears showing where are my "python paths" and after that I resolve put the script in "C:\OSGeo4W\bin" but I am not sure if here is the correct local.

The photo below show what happens when I try to load the shapefiles by python:

@gene, I follow your instruction, but how I update the question a mistake persists in the running process. First I have many "python paths" and I don't know where do I have to put the script, second when I use the command ldr.load_shapefiles('F:/qgis_1/exercicio_buffer_L1') I recive an error message like this in the photo.

What do I have to do?

///////////////////////////////////////////////////////////////////////////////////////////


texts:
I am running the following code and get:

Segmentation fault (core dumped)

I have no clue why, I have been using the functions many times before and now suddenly this happens. Sometimes it happens at the first read of ReadAsArray, and somtimes after 20+ reads.

#!/usr/bin/python3
import numpy as np

from methods import sampler

maps = [
    '{}_slope.tif'.format(continent) for continent in ['au', 'eu']
]

def create_tiles(size):
    for ulLat in range(90, -91+size, -size):
        for ulLon in range(-180, 181-size, size):
            yield ulLon, ulLat, ulLon + size, ulLat - size

def sample_point_xy(px, py, rb, cols, rows, win_size=1, func=None):
    if (px >= 0) & (px <= cols) & (py >= 0) & (py <= rows):  # check if within map extent
        if win_size > 1:
            intval = rb.ReadAsArray(
                int(px - win_size / 2.), int(py - win_size / 2.), win_size, win_size)
        else:
            intval = rb.ReadAsArray(
                int(px), int(py), win_size, win_size)
        if func is not None:
            value = func(intval)
        else:
            value = intval
    else:
        value = np.nan
    return value

def get_samples():
    size = 4
    for m in maps:
        rb, gt, cols, rows = sampler.get_raster_info(m)
        for tile in create_tiles(size):
            ulLon, ulLat, lrLon, lrLat  = tile
            ulx, uly = sampler.MapToPixel(ulLon, ulLat, gt)
            lrx, lry = sampler.MapToPixel(lrLon, lrLat, gt)
            win_size = int(max(lrx-ulx, uly-lry))
            px, py = (ulx + lrx) / 2, (uly + lry) / 2
            sample = sample_point_xy(px, py, rb, cols, rows, win_size=win_size)
            yield sample

for sample in get_samples():
    if sample is not None and not np.isnan(sample):
        print(sample)
        print(sum(sum(sample)))


texts: MapServer: Output angled label with Cairo/PDF

texts:
We're trying to output a PDF with the Cairo driver.
When printing a label with ANGLE set to follow, auto or an attribute, it fails regularly.

Has anyone experienced this before? A solution or workaround?

/Kristian


texts: Building a grid of equal square polygons in ArcGIS 10?

texts:
Not having them depend on decimal degrees, but on a certain distance.

When I use a tool it asks me the length of the cell in degrees, but I want all sides of every cell to be equal and with degrees it always makes rectangular cells or different lengths.


texts:
ogr2ogr -overwrite -t_srs EPSG:3978 -f "ESRI Shapefile" -sql "SELECT * FROM 'Layer Name' WHERE FIELD_NAME = 'A CLASS'"  dst.shp src.shp

help me out here - looking to
1) reproject a shape and
2) select only "A CLASS" from field named "FIELD_NAME"
3) and save it to src.shp

my output is
ERROR 1: Type mismatch or improper type of arguments to = operator.

I've also tried
ogr2ogr -overwrite -t_srs EPSG:3978 -f "ESRI Shapefile" -where "FIELD_NAME = 'A CLASS'"  dst.shp src.shp
resulting in:
ERROR 1: SQL Expression Parsing Error: syntax error
FAILURE: SetAttributeFilter( WHERE FIELD_NAME = 'A CLASSS') on layer 'LAYER NAME' failed.

and

ogr2ogr -overwrite -t_srs EPSG:3978 -f "ESRI Shapefile" -where "FIELD_NAME LIKE 'A CLA%'"  dst.shp src.shp
with the las
which only produces one results, where 10,000+ records satisfy the query.

update

ogr2ogr -overwrite -t_srs EPSG:3978 -f "ESRI Shapefile" -dialect SQLite -sql "SELECT * FROM 'Layer Name' WHERE FIELD_NAME = 'A CLASS'"  dst.shp src.shp

yields
ERROR 1: Type mismatch or improper type of arguments to = operator.
ERROR 1: In ExecuteSQL(): sqlite3_step(SELECT * FROM 'Layer Name' WHERE GIVEN_CLASS = 'Class A'):
  SQL logic error or missing database

one step closer though!


texts: ArcGIS - Show Point Marker Symbol based on attribute

texts: I am using ArcGIS 10.4 and I have a point layer of universities. I want to show their logos in map layout. For each record in attribute table I have the link to respective logo of that university. I am using an option to assign a point marker symbol but the issue is it is only allowing me to pick one by one by browsing the file explorer and choosing an image for each point. I am not getting a way so that it picks from attribute table and then show the image.


texts: Weighted Overlay - passing parameter as text error?

texts:
For some reason when hardcoding the save path for the wieghted raster just before runtime it works perfectly.  Like such...

arcpy.env.cellSize = "90"
arcpy.env.workspace = r'G:\Xcel\Route Tool\Potter\Scratch.gdb'
landuseRaster = r'G:\Xcel\Route Tool\Potter\Data.gdb\Potter_Randall_Landuse_Raster'
slope_Reclass_Raster = r'G:\Xcel\Route Tool\Potter\Scratch.gdb\slope_Reclass'

#Generate a list identifying what the individual input values should be reclassified as.
remapLanduse = RemapValue([['Cropland and pasture',2],['Herbaceous Rangeland',1],['Nonforested wetland',10],['Other agricultural land',1],['Strip mines, quarries and gravel pits',2],
                            [ 'Mixed rangeland',1],['Lakes',10],['Residential',3],['Commercial and Services',3],['Confined feeding operations',3],['Transportation, communications and services',10],
                            ['Transitional areas',3],['Shrub-brushland rangeland',1],['Other urban or built-up land',3],['Bare exposed rock',3],['Industrial',2],['Reservoirs',10],['Forested wetland',10],['Mixed urban or built-up land',3]])

remapSlope = RemapValue([[1,1],[2,2],[3,3],[4,4],[5,5],[6,6],[7,7],[8,8],[9,9],[10,10]])

#Define the input rasters, the field identifying the input values, the remap of their values, the weight of each raster, and the evaluation scale to use in the Weighted Overlay tool.
myWOTable = WOTable([[landuseRaster, 80, "LANDUSE", remapLanduse],
                     [slope_Reclass_Raster, 20, "Value", remapSlope],
                     ], [1, 10, 1])

outWeightedOverlay = WeightedOverlay(myWOTable)

arcpy.AddMessage("Creating weighted raster...")
weighted_Raster = env.workspace + '\\weighted_Raster'
outWeightedOverlay.save(weighted_Raster)

However when I try to pass the workspace to weighted overlay tool using parameter as text I receive a generic 999999 error? Like such....

subs  = arcpy.GetParameterAsText(0)
Field = arcpy.GetParameterAsText(1)
originValue = arcpy.GetParameterAsText(2)
destValue   = arcpy.GetParameterAsText(3)
lines = arcpy.GetParameterAsText(4)
Dem = arcpy.GetParameterAsText(5)
landuseRaster = arcpy.GetParameterAsText(6)
outwork = arcpy.GetParameterAsText(7)

arcpy.env.workspace = outwork + "\\Scratch.gdb"

.... run a few tools including generate 'slope_Reclass_Raster'

#Generate a list identifying what the individual input values should be reclassified as.
remapLanduse = RemapValue([['Cropland and pasture',2],['Herbaceous Rangeland',1],['Nonforested wetland',10],['Other agricultural land',1],['Strip mines, quarries and gravel pits',2],
                            [ 'Mixed rangeland',1],['Lakes',10],['Residential',3],['Commercial and Services',3],['Confined feeding operations',3],['Transportation, communications and services',10],
                            ['Transitional areas',3],['Shrub-brushland rangeland',1],['Other urban or built-up land',3],['Bare exposed rock',3],['Industrial',2],['Reservoirs',10],['Forested wetland',10],['Mixed urban or built-up land',3]])

remapSlope = RemapValue([[1,1],[2,2],[3,3],[4,4],[5,5],[6,6],[7,7],[8,8],[9,9],[10,10]])

#Define the input rasters, the field identifying the input values, the remap of their values, the weight of each raster, and the evaluation scale to use in the Weighted Overlay tool.
myWOTable = WOTable([[landuseRaster, 80, "LANDUSE", remapLanduse],
                     [slope_Reclass_Raster, 20, "Value", remapSlope],
                     ], [1, 10, 1])

outWeightedOverlay = WeightedOverlay(myWOTable)

arcpy.AddMessage("Creating weighted raster...")
weighted_Raster = env.workspace + '\\weighted_Raster'
outWeightedOverlay.save(weighted_Raster)


texts: Multivalue iteration in Vehicle Routing Problem (VRP)?

texts: I'm currently building a model that would simulate three different set of orders (shapefile format) by using the Vehicle Routing Problem (VRP) tool of network analyst. The model I built runs correctly but shows me only the first delivery solution in the excel file. My aim is to have all the different route results into an Excel file (or even multiple). Does anybody know how to do that or fix the model? I also add that one issue is related to the fact that I cannot specify the input data in the delivery "Add location" during different iterations (only the first one is considered).


texts: Can't connect to WMS service? An Amazon-EC2 problem?

texts: I'm running GeoServer 2.3.5 on a Windows instance in AWS for the first time.  The data is on the server in PostGIS.  The map layers have been published and I can see them locally without any problems.  However, when I try to access the WMS service for these layers from another machine, I can't connect.  I'm wondering if there is some configuration I may be missing because this is set up through AWS.  Any ideas?


texts: Automatically replace values in Geometry generator's expression in all symbol layers of feature

texts:
I have created some piecharts in QGIS 3.14 with the geometry generator (see here why it was not possible to use the built-in piecharts from the layer properties).
I now have 10 piecharts, each with several sectors (slices) styled in different pattern-styles, see screenshot below for how it looks. Some of the piecharts are not located in an ideal place, so I want to re-position them somewhere else - using the project(point,distance,azimuth) expression in the geometry generator.
However, for this I have to change the expression for each sector/slice separately as I built the piechart with wedge_buffer, creating a separate symbol layer for each sector/slice.
I could use a variable as explained here, but that's not the point since I would have to introduce the variable manually for each symbol layer.
Is it possible to automatically replace the values of distance and azimuth (or any other part of the expression in the geometry generator) automatically, by iterating over all symbol layers (like the highlighted examples no. 1 and no. 2 in the screenshot) and not by changing each symbol layer manually one by one? I am looking for something like "search for x and replace with y" over all symbol layers of a feature.


texts:
I am looking for a way to transform a QGIS project (.qgs) to a MapServer Mapfile (.map).

I am using Quantum GIS 2.6 and MapServer 6.

The layers are postgis tables.

I would need something like this but it is not being displayed in my QGIS plugins list: http://docs.qgis.org/1.8/es/docs/user_manual/plugins/plugins_mapserver_export.html

Instead it only shows this plugin: https://plugins.qgis.org/plugins/rt_mapserver_exporter/ but does not seem to convert qgs to mapfile.

Any solution to this?


texts:
In QGIS we have new layout type: Report. I didn't find tutoral to export this in PDF in Python.
Where can I find Python code examples?


texts: QGis [2.14 / 2.18.10] crashes selecting a feature [sqlite database]

texts:
I have a QGis Project with five vector-layers; every layer has about 100 up to 400 features. Everything is saved in a sqlite-database.

When I take the selecting tool (Select Features by area or single click), QGis crashes in the moment when I select a single item. Sometimes the same happens when using the Node Tool in the same moment. I can't recognize a logic to reproduce this failure. I can say that it happens not so often when saving the layer and the project frequently (and switch the Toggle Editing on and off). But it appears always the same error message:

  qgis.bin: geos_ts_c.cpp:3472: int
  GEOSCoordSeq_setOrdinate_r(GEOSContextHandle_t,
  geos::geom::CoordinateSequence*, unsigned int, unsigned int, double):
  Assertion `0 != cs' failed.

Tested in QGis 2.14 and 2.18.09 / 2.18.10


texts: Georeferencing and mosaicking jpeg images using ArcMap with unknown spatial reference?

texts:
I have a series of large scale (almost 1:1), aerial jpeg images taken from a UAV DJI-FC350 drone camera. They have no spatial reference, but the images have their own GPS tags giving a single lat, lon and altitude for each.

I am trying to see if it is possible to use ArcMap to load them as raster datasets and mosaic them together to create a single image with the proper spatial locality.

I have read in some places that this might not be possible and you must use a photo stitching program to get a single image.


texts: Can AutoCAD connect to PostGIS (specifically for Ordnance Survey Mastermap)?

texts:
As per title.

Can up to date (2019) AutoCAD versions be practically/safely connected to a PostGIS database.

And more specifically does this work well/practically/safely for accessing Ordnance Survey Mastermap data? What are the limitiations?


texts: Calculating distance between polygons using ArcMap

texts:
I am a new user to ArcMap.
I have two layers of polygon features and I want to calculate the minimum and mean distance between polygons of one layer and polygons of the other layer.
I am guessing that Euclidean Distance is a tool I would use, but I'm not positive.
How do I approach this problem?


texts:
I'm using Agisoft Photoscan 1.4.2 to build DEM from UAV photos and in last versions thay produce DEM multipage tiff format and I want to convert it to tiff with metainformation preserved.

According to this, I use vips to convert multipage tiff to tiff, but it drops all metainformation:

sudo apt-get install libvips-tools
vips im_vips2tiff dem_1_ctb.tif:0 1.tif

Then I wrote python script using GDAL api:

from osgeo import gdal, gdalconst

src = gdal.Open('dem_1_ctb.tif', gdalconst.GA_ReadOnly)
dst = gdal.Open('1.tif', gdalconst.GA_Update)

prj=src.GetProjection()
print('prj:', prj) #

gtr = src.GetGeoTransform()
print('gtr:', gtr) #

dst.SetProjection(prj)
dst.SetGeoTransform(gtr)

Original metainformation obtained by gdalinfo:

Driver: GTiff/GeoTIFF
Files: dem_1_ctb.tif
Size is 1773, 1603
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (31.826655569140200,59.242815662671902)
Pixel Size = (0.000002169874100,-0.000001111660512)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  31.8266556,  59.2428157) ( 31d49'35.96"E, 59d14'34.14"N)
Lower Left  (  31.8266556,  59.2410337) ( 31d49'35.96"E, 59d14'27.72"N)
Upper Right (  31.8305028,  59.2428157) ( 31d49'49.81"E, 59d14'34.14"N)
Lower Right (  31.8305028,  59.2410337) ( 31d49'49.81"E, 59d14'27.72"N)
Center      (  31.8285792,  59.2419247) ( 31d49'42.88"E, 59d14'30.93"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  NoData Value=0
  Overviews: 887x802, 444x401, 222x201
  Unit Type: metre

Transfered metainformation:

Driver: GTiff/GeoTIFF
Files: 1.tif
Size is 1773, 1603
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (31.826655569140200,59.242815662671902)
Pixel Size = (0.000002169874100,-0.000001111660512)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=3 (pixels/cm)
  TIFFTAG_XRESOLUTION=10
  TIFFTAG_YRESOLUTION=10
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  31.8266556,  59.2428157) ( 31d49'35.96"E, 59d14'34.14"N)
Lower Left  (  31.8266556,  59.2410337) ( 31d49'35.96"E, 59d14'27.72"N)
Upper Right (  31.8305028,  59.2428157) ( 31d49'49.81"E, 59d14'34.14"N)
Lower Right (  31.8305028,  59.2410337) ( 31d49'49.81"E, 59d14'27.72"N)
Center      (  31.8285792,  59.2419247) ( 31d49'42.88"E, 59d14'30.93"N)
Band 1 Block=1773x16 Type=Float32, ColorInterp=Gray

As I can see some of the metainformation is different, is it harmful? can it be transferred completely?


texts: How to change the dynamic information in data driven pages

texts:
I want to ask one question about data driven pages.

I plan to create data driven pages to show each parish’s information in Louisiana. Followed is the template I am going to use. I downloaded it from ATSDR.
See there is Demographic Statistics in the plate? I wish the statistics could change when the county changes.
For example, when the map demonstrates county A, county A’s information would appear; while in next page, the map showed county B, and county B’s information come out.
I tried to insert dynamic elements, but it failed. Because dynamic elements only include page number, date, etc.
Is there anyone know how to deal with that? If you need further information, please let me know.


texts: Building an interactive web map that displays table data in a static portion of the page separate from the map?

texts:
I am trying to build an interactive web map to display the results of my master's thesis. Basically, it is a shapefile of line segments that have each been assigned to 1 out of 8 categories. The category is contained in the shapefile's attribute table.

When a user clicks on a particular line segment, I would like a static part of the page, separate from the map, to be populated with information based on what category the line segment is in.

I know it is fairly easy to set up a map in MapBox or Carto where you can click on segments and a little bubble pops up with the data contained in the attribute table. However, the information I want to display would be in the form of a rather large HTML table (1 for each category) that would best be displayed in static portion on the page to the side of the map, rather than a pop up bubble.

I am familiar with Mapbox, Carto, QGIS2Web, and am currently learning HMTL/CSS and JavaScript. Right now I am planning to build the basemap in Mapbox, export the shapefile as a geojson Leaflet map using QGIS2Web, and then code the functionality I want for the static portion of the page in Javacript. The HTML code for the table would be far too large to store in the attribute table of the shapefile, so I would have to write a script that would display one of the 8 tables based on the category code for a given segment.

But before I get too far into the weeds, does anybody know of a mapping platform that allows for this functionality? Or a way to do it that requires less custom coding?

I'm happy to put the work in, but would kick myself if I found out there was an easier way to do it afterwords.


texts: Determining which specific features are visible from observer points using ArcGIS 3D Analyst?

texts:
I am working on a wind farm proposal. I have the positions of 70 proposed turbines. I also have the positions of 80 houses within 20km of the proposed turbines. I created a viewshed for the 70 turbines, and was able to determine which parts of the surrounding landscape would be able to see at least one turbine. The viewshed has also told me how many turbines are visible from each pixel in the raster. Based on the position of the house, I have determined how many turbines are visible from each house.

But now I'd like to take it a step further. I would like to determine which specific turbines are visible from each house.

I have two ideas, both of which involve a lot of labor. (1) Create a viewshed for each house individually. Record the specific turbines visible from each house. This will involve creating 80 individual viewsheds - too many. (2) Use the "Create line of site tool" between each house and each turbine. This will involve 80 x 70 = 5600 individual calculations - too many.

There must be a quicker way to do this. I have ArcGIS 10.0 with 3D analyst activated.


texts: Zonal Statistics of a polygon and assigning mean value to the polygon

texts:
I have a large shapefile filled with polygons, and a raster of that same area. I need to use Zonal statistics to find the mean value of the area under each polygon and then assign that value to the polygon.

What would be the best way to assign the mean values to their corresponding polygons for later retrieval?

I am working in Python, and using mainly GDAL/OGR.


texts:
I am still relatively new to QGIS and have been using it for the basemap image services associated with OpenLayers plugin_particularly, the Google Satellite service. I am now unable to add/see the All plugins option in the QGIS GUI.

I read the problems in the context of OL so I'm trying to access QuickMapServices but I'm not able to see any of the plugin options except those which have been installed in my previous version_now updated to Essen. Oddly, and pertinent to this question I think, is the fact that I downloaded the software on a new laptop and was able to access all plugins and see them successfully.


texts: How do I add elevation data to a GeoJSON line from a GPX file?

texts: I understand that elevation can be added as an optional third value in Coordinates (lat, long, elev) in a geojson file, but I can't find a way to actually do that from a gpx file. I've used ogr2ogr to convert gpx to geojson.  I've also loaded the gpx file into QGIS and saved as a geojson file, but neither method brings in elevation when I do this to the gpx tracks. It will bring in elevation as a Property if I convert the points, but I need a line. Ultimately, I want to diplay these lines in Leaflet and build elevation profiles for them.  Is there any easy way to do this other than parsing out the gpx files and building my own geojson files?


texts: Python syntax error...cannot see the error!

texts:
I'm very new to python and have created a script for enabling a user to list file types in a workspace making use of the parameters function in ArcGIS so it's more user friendly. However I cannot get the script to work...I am struggling to finish the script off can anyone lend a hand?

Python reports the syntax error in line 36 "f.close()" there is no white space around it so am struggling to see whats wrong...can you? Or what am I missing?

Here's the code:

import arcpy
from arcpy import env
import os
env.workspace = ("C:\ArcGIS\Default.gdb")
ws = arcpy.GetParameter(0)
option = arcpy.GetParameterAsText(1)
outputFilePath = arcpy.GetParameterAsText(2)

dsList = []

if option == "List Feature Classes":
    featureclasses = arcpy.ListFeatureClasses(ws)
    for table in tables:
        dsList.append(table)
elif option == "List Tables":
    tables = arcpy.ListTables(ws)
    for table in tables:
        dsList.append(table)
elif option == "List Rasters":
    rasters = arcpy.ListRasters(ws)
    for table in tables:
        dsList.append(table)
elif option == "List Workspaces":
    workspaces = arcpy.ListWorkspaces(ws)
    for table in tables:
        dsList.append(table)
elif option == "List Files":
    files = arcpy.ListFiles(ws)
    for table in tables:
        dsList.append(table)

f = open(outputFilePath, 'w')

for item in dsList:

f.close()


texts: Scale and Z factor have no effect on hillshade analysis in QGIS

texts:
I used an SRTM DEM to create a hillshade in QGIS 2.2.0. However, the resultant hillshade seems to be vertically exaggerated. Also, changing the settings for scale and z-factor does not have any effect on the resultant hillshade. I have tried both options in QGIS, Analysis-> DEM (Terrain models) and Terrain analysis -> Hillshade. Both methods gave the same results.

Is this an isolated incident, a bug, or am I missing something?


texts: A plugin or script for separating the elements of a raster

texts: I am looking for a plugin or script that allows  to separate the elements of a raster (photovoltaic panels in a solar power plant), each in an image.


texts: Can I get a .jpg from a CartoDB map created in the editor without using a browser?

texts:
Edit

I'm actually starting to suspect that I can't get a .jpg image from a map created with the UI editor without using a browser.  I think I might have to use the API to re-create the map my colleague created so that the map has a layergroupid.  An acceptable answer to this question might be to confirm those suspicions.

Original question

I'm a complete beginner with CartoDB.

Somebody in our organisation has created a map in the CartoDB editor (described here), and I need to get a .jpg image from it.  The map covers the whole of the country and shows property boundaries.  The person who created it is responsible for the styling of the map and for the data the map uses.  The datasets are private.

I am writing an API which needs to get aerial images of properties from the CartoDB map, zoomed in so that each image shows a few properties.  My API is written in C# and runs on a server.  I mention those details to emphasise that there is no browser involved in the interactions with the map, and no Javascript.

The map need only be accessible to trusted entities.  The two I can think of are: (1) the person maintaining the map and (2) my code.  The code I'll be writing runs in a totally private context on our server.  It would therefore be OK for my code to use the api_key when retrieving the .jpg.  I can also set a password for the map and use that.  (The images retrieved will be stored on our company's server, and those stored images served up to our API users in a completely different security context.)  So access to the map is totally private.

I've seen some information about getting .jpg images, but it all seems not to work for me because it requires a layergroupid, which it seems is only available if the map was created through a CartoDB API.  Since my map was created using the editor, I don't think there is a layergroupid, and so those methods don't apply.

It might be possible to create the map using the CartoDB API instead of using the editor, so that the map has a layergroupid.  This approach has drawbacks though, because my colleague would need to learn a different way to create the maps.

In summary

complete beginner with CartoDB
map created using CartoDB editor, private datasets
no browser, no Javascript
method of retrieving .jpg totally private, so could use api_key
prefer not to re-create existing editor-created map via CartoDB API

Question

Given those constraints, how can I get the .jpg images I need from the map?


texts: Define vector projection and reproject shapefile in ERDAS IMAGINE 2013?

texts:
I come from a GIS background mainly using ArcMap and QGIS, but am now doing more image processing work and so have decided to start using ERDAS IMAGINE 2013.

I am attempting to project a shapefile from WGS 1984 to British National Grid.

I open the shapefile, and it displays fine:

From the vector tab, I select "Reproject shapefile". For the input file I select this shapefile I opened, however I get the error "Fail to convert string to coordsys obj." Googling this error produces no results. The source projection fields are not filled, and I cannot reproject the shapefile:

Examining the metadata for the shapefile, I notice the projection tab is blank:

Selecting Edit -> Add coverage projections (I assumed this was a way of defining the projection for the layer), I select Edit Projection Parameters, and then select WGS 84 and Decimal Degrees. I select Yes overwrite PRJ file, then get the error message "Faled to translate erdas map projection to PE coordsys.", which again brings up no results in Google.

ERDAS seems to be unable to read projection information for my shapefiles, or possibly unable to read the .prj files. Watching this video using an older version of ERDAS: https://www.youtube.com/watch?v=Vfk4eOJ92xw the process seems straightforward enough, but I've tried various shapefiles from different sources, and new shapefiles made in Arc and QGIS, and get the same error. I've managed to reproject raster data fine.

I can reproject the shapefiles fine in ArcMap, but like I said I'm trying to get used to ERDAS as much as I can.

Is there a step I'm missing somewhere?


texts: Problem exporting a large classified image in Earth Engine

texts:
I am trying to export a GeoTIFF image but it doesn't work. https://code.earthengine.google.com/fe534ace7f3fc4fdfd155ded24155525

When I export the large classified image to Google Drive, I don't get an error and the task takes a long time and finally, I can't see the GeoTIFF image.

Here is the export request:

 Export.image.toDrive(
                {image:PBS_OUT.float(),    //'B2','B3','B4',
                description: 'PBS_S2_2019',
                fileNamePrefix:'PBS_S2_2019',
                 folder: '',
                 region:AOI.geometry().bounds(),
                 scale:20, //10
                 skipEmptyTiles:true,
                 maxPixels: 1e13
                });

Any ideas for getting my export task to be successful?


texts: Filter layers on layer group with OpenLayers and GeoServer?

texts:
I have a layer group with multiple layers as following

In Angular I created an OpenLayers map adding the group layer like this:

this.capaPropia =
new ImageLayer({
  source: new ImageWMS({
    url: environment.baseUrl + 'handlerGeoserver.ashx',
    params: {LAYERS: 'v_expedientesAgrupados'},
    projection: proyeccion
  })
});

In my app I have 4 checkboxes, one for each layer. I need to display only the layers of the group layer that are checked.

Is it possible to filter layers from a group layer on GeoServer (with cql_filter or similar)? If so, can this be done with OpenLayers?


texts: PostGIS 2.1 ST_EstimatedExtent returns error

texts:
In absence of statistics (empty table or no ANALYZE called) the function ST_EstimatedExtent(schema,table,geom_column) should return NULL (PostGIS 2.1).

In my case, the function returns

FEHLER:  stats for "table1.the_geom" do not exist

********** Fehler **********

FEHLER: stats for "table1.the_geom" do not exist
SQL Status:XX000

The table table1 is empty, so the function should return NULL?!
I also run an ANALYZE on the table.


texts: I have an attribute field which contains Date data. The data is broken down like "3/14/2016". I would like to have different symbology for different years. When I go into the symbology tab within the layer properties and set the field that contains the date I do not have the option to set the symbology for each year. Instead it shows me a list of dates contained in the fields like the above example. How do I tell arcmap to only look at the year and not the month and day? So in the long run every feature with the same year has the same symbol?


texts: How do I change the color of an ArcGISImageServiceLayer?

texts:
I have published an image service of a raster layer (.tif).
When I use it online, it is grayscale.
In ArcMap, I have the option to change the color of the display.
How can I do this from the webservice (either when publishing, or consuming)?

Here is the code I use to consume the service.

var params = new esri.layers.ImageServiceParameters();
    params.noData = 0;

var montgomeryCountyLayer = new    esri.layers.ArcGISImageServiceLayer(".../arcgis/rest/services/Montgomery_Layer/ImageServer", {
      imageServiceParameters: params,
      opacity: .5,
      id: "montgomeryCounty"
    });

    map.addLayer(montgomeryCountyLayer);


texts: Getting screen coordinates from canvas coordinate using PyQGIS

texts:
I'm trying to set the initial position of my QGIS plugin dialog. I'm using .setGeometry() but the problem is that I don't know how to get coordinates. I want to set the initial position in the top corner of the canvas in QGIS (like the image below). I've tried this but it does not work.

self.dlg = SetZoomDialog()
canvasCenter = iface.mapCanvas().extent().center()
canvasToGlobal = iface.mapCanvas().mapToGlobal(QPoint(int(canvasCenter.x())/2, int(canvasCenter.y())/2))
x= canvasToGlobal.x()
y= canvasToGlobal.y()
self.dlg.setGeometry(x,y,width,height)

I tried also this canvasToGlobal = QgsMapTool(iface.mapCanvas()).toCanvasCoordinates(canvasCenter) but still I don't understand how to convert canvas coordinates to screen coordinates.


texts: Writing expression for field calculator in Python?

texts:
I have a script where I calculate a couple of variables, add a new field then need to do a field calculation that field. Once I get everything set up, it's the expression that I have a hard time figuring out. I am using ArcGIS 10.2

denominator = 0.14
arcpy.CalculateField_management("My_Table", "LQ", ( ( "!Foreign_Born_Population!" / "!Population!" ) / denominator ) , "PYTHON_9.3")

I am thrown an error stating it cannot parse the string but I thought by wrapping the field name in !! that it would recognize it as a field.

If I were to open the field calculator in the attribute table and type in a function using VB Script, the function would look like:

( [Foreign_Born_Population] / [Population] ) / 0.14

I also tried using the VB Script in the function but that did not work.

I have seen a couple different solutions including with and without using a cursor. Do I need to use a cursor?

How can I write this function?


texts: Options for ArcGIS Server and ArcSDE Hosting

texts: I know this was partially address here, but it was not specific to ESRI products.  I want to find some place to create and host map services with ArcGIS Server and ArcSDE.  ESRI cloud options are bit pricey for the initial startup and testing of my project.  Anyone know of lower priced alternative?


texts: Mapping objects/assets along polygon line using ArcGIS Pro?

texts:
I am working with an open rail data set which is imported to my ArcGIS Pro project. The rail data is loaded to ArcGIS as a shapefile and includes kilometer notation for a specific rail segment. The km notation is included for every 0.1 km.

I want to add new object/assets but I only have the km notation as an input and need to align them along the curve rail track as visually shown below.

Could someone guide me how to map the new objects?


texts: Cannot establish connection to port 6080

texts:
From what I understand, the default port for ArcGIS for Server is port 6080. I have installed ArcGIS for Server 10.2 and now I cannot connect to the Server Manager. I am using Windows 8. My firewall is completely turned off. I have restarted the service about 100 times and I have reinstalled ArcGIS for Server about 10 times at this point. The port is not shown when I do a netstat -an in command prompt as well.

URL I used: http://localhost:6080/arcgis/manager

Picture of "error":


texts: Can ArcMap basemaps be used in figures intended for publication?

texts: I'm wondering if the 'imagery' basemap available in ArcMap 10 can be used as part of a figure in a scientific publication? If it can't are there any alternatives, using QGIS for example?


texts: Is there any method to create monthly LST images?

texts:
Is there a method to create monthly LST images from combining MOD11A2(8days) product?

The MOD11C3 product which is delivered monthly has larger spatial resolution than MOD11A2.

Im trying to create a script in python to do that but it seems to be above my coding skills


texts: Separating polygons within layer and measure line through each polygon using ArcGIS Desktop?

texts:
I have 3 layers of polygons, each layer with multiple polygons spatially separated, although the layer only counts it as 1 polygon.
I have a layer with a line going through these polygons.

What is the best way to measure the line running through each polygon? I would like a measurement for each individual polygon in each layer, even though each layer is only considered as one whole polygon, I want it to be separated. I prefer not to do it manually, as there are over 500 polygons.

Above is a picture of a section of the map I have so far.

I'm using ArcMap 10.1.


texts:
I use ArcGIS 10.1 for Desktop.

I have a polygon Feature Class with municipalities.

In a non-geographic table, I have these same municipalities, with an additional contract information.

The contract name is stored as a String in a field. Each municipality can have either zero, one or more contracts. In the table it appears as often as it has the information, i.e. a municipality with three contracts appears three times in my table.

At the end, I need a dissolved polygon for each contract where every munipality of the specific contract is within. These polygons must be overlapping when a municipality has more than one contract.

Does anyone have any idea how to to this, preferably in ModelBuilder, that another user could reuse it, only pointing where he has the existing table and feature class?


texts: QGIS update feature geometry from attribute fields

texts:
I know how to get the geometry (x,y) into attribute fields with the calculator, but I can't figure out how to update the geometry from values contained in attributes.

I am trying Field Calculator > Update existing field > , but I can't find how to assign my "x" and "y" fields to the actual feature geometry.

I don't want to manually drag each point with the vertex editor, or type in each coordinate. I have to update many points.


texts: Compute neighborhood statistics in GEE with your own function/reducer

texts:
I have generated an Image in GEE, which corresponds to a classification. Each pixel is assigned to 1 group (out of a certain number of groups, let’s say 100). I would like, for each pixel, to compute neighborhood statistics within a focal window of a given size, say a square window of size 10 pixels. I found the “reduceNeighborhood” function which does just that. However, there is only a handful of reducers available, and the neighborhood statistics I’m trying to compute are slightly more complex than what the available reducers allow doing. For instance, I’d like to compute a diversity index, such as the Shannon index. This is a very simple formula (see Wiki), but I cannot compute it by combining existing reducers (at least I hardly see how).
Given my limited knowledge of GEE, I foresee two approaches to compute the Shannon index:
1/ The first one is the simplest one: specify my own formula to be used as a reducer in the “reduceNeighborhood” function. I’ve looked around but never saw someone do that, so I’m guessing it’s not possible.
2/ The second one is the painful one: (1) disaggregated the classified image into 100 images (1 per group), (2) compute some neighborhood statistics on each of the 100 images (i.e. the abundance – or count – of pixels belonging to each group within focal windows), (3) reaggregate the 100 output bands into 1 using the appropriate formula (Shannon index). This strategy would require a loop and seems like a very inefficient use of computational resources.
Thus, before trying to figure out if 2/ is implementable, I’d like to make sure option n°1 is indeed not possible. My question this is : Is it possible to specify its own function / reducer in “reduceNeighborhood” ?


texts: Python Script to apply the Northing and Easting of a photo center .shp file to a tiff.

texts: I have a group of tiffs I want to georeference. I have the photo centers for each of the tiffs but they are not linked. After looking at the comments I believe there would be a way to take the attribute of the .shp, Northing and Easting, and apply it to the tiff? The frame number for the photo is in the attribute table of the points. I assume you could use this to look into the folder and pull out the correct photo and reference it to said location. I hope this explains what I am looking for in a better light.


texts: Overpass XML query for finding all buildings inside a polygon

texts:
I want to run following Overpass QL query in order to get way IDs of all buildings inside a polyon:

way(poly:"50.7 7.1 50.7 7.2 50.75 7.15")[building];
>;
out ids;

Since I'm using the osm-common library, I need to convert this script to Overpass XML. I did it using the convert form here.

The result:

<osm-script>
    <query into="_" type="way">
        <polygon-query bounds="50.7 7.1 50.7 7.2 50.75 7.15" into="_"/>
        <has-kv k="building" modv="" v=""/>
    </query>
    <recurse from="_" into="_" type="down"/>
    <print from="_" limit="" mode="ids_only" order="id"/>
</osm-script>

But when I run it, I get node IDs, not way IDs.

<?xml version="1.0" encoding="UTF-8"?>
<osm version="0.6" generator="Overpass API">
<note>The data included in this document is from www.openstreetmap.org. The data is made available under ODbL.</note>
<meta osm_base="2013-11-06T07:51:02Z"/>

  <node id="252553109"/>
  <node id="252553111"/>
  <node id="252553112"/>
  <node id="252553114"/>
  <node id="252553115"/>
  <node id="252553117"/>
  <node id="255302992"/>
  <node id="255302993"/>
  <node id="255302995"/>
  <node id="255763355"/>
  <node id="255763356"/>
  <node id="255763357"/>
  <node id="255763358"/>
  <node id="255763430"/>
  <node id="255763431"/>
  <node id="255763432"/>
...
  <node id="2519505892"/>
  <node id="2519505893"/>
  <node id="2519505894"/>
  <node id="2519505895"/>
  <node id="2519505896"/>

</osm>

How should I modify the Overpass XML query in order to get way IDs of the buildings?


texts: How configure ArcGIS 10 Server for .NET at http://server/rootContext/ArcGIS rather than http://server/ArcGIS

texts:
We are developing a Flex application which will communicate with ArcGIS server via a reverse proxy. The reverse proxy configuration is such that the host and path portions are separate configurations. Does anyone know if you can configure ArcGIS to respond to a URL like:

http://serverName/rootContext/ArcGIS

rather than the default:

http://serverName/ArcGIS

Has anyone done this? Any suggestions on how that might be done?

Thanks,
Carl


texts: Wrong coordinates of mapproxy TMS in Leaflet 0.7

texts:
My problem is similar to this question. The coordinates of my Mapproxy cache are wrong: Greenwich has latitude 180 instead of 0. See the picture below. The option tms: true was used as you may see from the code, but the boundaries of the map are still wrong. What should I do to have correct coordinates of the cache?

import L from 'leaflet';

class MapController {
constructor(Map, $window, intercom) {
    'ngInject';

this.$window = $window;
this.Map = Map;
this.intercom = intercom;

this.Map.status = "opened";

// Закрытие окна карты при закрытии окна приложения
this.$window.onunload = () => {
this.Map.status = null;
};

this.intercom.on('map-event', this.mapEventHandler.bind(this));

var osm = L.tileLayer('http://gdb-dev:8080/mapproxy/tms/1.0.0/OSM_4326/EPSG4326/{z}/{x}/{y}.png', {tms: true});
var bing = L.tileLayer('http://gdb-dev:8080/mapproxy/tms/1.0.0/Bing/EPSG4326/{z}/{x}/{y}.png', {tms: true,continuousWorld: true});
var googlemercator = L.tileLayer('http://bankrl:8080/mapproxy/tms/1.0.0/Google-mercator/EPSG900913/{z}/{x}/{y}.png', {tms: true });

let map = new L.map('map', {
  layers: [osm, bing, googlemercator]
});
this.map = map;

map.on('load', () => {
console.log('Map is ready!');
this.Map.status = 'ready';
this.intercom.emit('map-ready', {message: "Hello!"});
});

osm.addTo(map);

var baseMaps = {
    "OSM" : osm,
    "BING" : bing,
    "Google-mercator" : googlemercator
};

L.control.layers(baseMaps).addTo(map);

L.control.mousePosition().addTo(map);

map.fitBounds([
    [-180.0, -90.0],
    [180.0, 90.0]
]);

map.setView(map.getBounds().getCenter(), 1);

}

}

export default MapController;


texts: WFS request using layer extent or spatial query QGIS3

texts:
I'm trying to request features from wfs-service by using a virtual layer. Is it possible to load features by using a layer extent from a virtual layer?

I know how to add dimension for a request with ogr2ogr but how can I do same thing in QGIS?


texts:
I am trying to work out what mapping tools I would need to use to dynamically create Choropleth maps from our data.

I am working on a project where the back-end is a big database of statistics, and the front end is various dynamically-generated presentations, from graphs to maps.

For example: If we had population statistics, we would want to be able to display them a little like this: http://crschmidt.net/mapping/choropleth.html or this: http://indiemaps.com/blog/2012/04/introducing-openlayers-symbology/ (although the examples for that library aren't working in my browsers).

Those are entirely in-browser: they take a GeoJSON file and adds it as a layer to OpenLayers. But it looks doubtful how responsive that can be: you seem to end up with either a fairly simple map which doesn't withstand zooming to country level well or a very large GeoJSON file. So I wonder if it is necessary to have some kind of server (MapServer?) dynamically generating map views at the requisite level of detail?

The locations in a dataset could be every country in the world, every country in a continent (or club, like OECD), or sub-national administrative areas. Ideally, as a version 2 aim, we want to be able to zoom in and out from world to sub-national level, as the project is all about providing more context for data.

EDIT: the Thematic Mapping API (http://thematicmapping.org/api/choropleth.php) looks like a solution but requires a browser plugin or KML viewer. We are looking for an in-browser solution.

So: Is dynamically generating Choropleth maps a solved problem? Is there a standard approach or library? If not, can you suggest a sensible angle of attack? I'm new to GIS and still reading in.


texts: How to display min and max pixel from Zonal Statistics?

texts:
I am trying to find the maximum and minimum pixel from a raster inside the area of a buffer. I can succesfully get the values using Zonal Statistics, but it doesnt tell me which pixel is the max and min and I would like to determine that in order to be able to display it in the map.

I am using QGIS 2.18.16


texts: Unable to use fillColor values from GeoJSON with Leaflet

texts:
I have a web page which fetches some remote GeoJSON and displays a combination of marker points, polygons and lines on my map.

Within that JSON is the fillColor, fillOpacity, etc that I wish to use for each polygon.

Here's an excerpt of the GeoJSON which is returned:

{
    "geometry": {
    "coordinates": [
        [
        [
            -2.4710206012402689,
            53.470236995186752
        ],
        [
            -2.4691196881249033,
            53.4613205079861
        ],
        [
            -2.4681978304127927,
            53.461390191450619
        ],
        [
            -2.4681448784765276,
            53.461140769261917
        ],
        [
            -2.4710206012402689,
            53.470236995186752
        ]
        ]
    ],
    "type": "Polygon"
    },
    "id": "CE71BB87D309C9199DC60BC7586774468D2AE99C",
    "properties": {
    "hazardFactor": "60",
    "hazardFactorName": "Danger",
    "fillColor": "#FF0055",
    "strokeColor": "#BF0442",
    "fillOpacity": "0.35",
    "strokeWidth": "1",
    "strokeOpacity": "0.8",
    "detailedCategory": "uk:frz",
    "iconUrl": "https://domain.com//v1/map/icon?icon=uk_frz.png",
    "name": "Manchester/barton",
    "pointOfContact": "0161-789 1362 (Administration and Operations)",
    "listOrderHint": "1000",
    "category": "airspace",
    "designator": "EGCB",
    "filters": [
        {
        "name": "Airspace",
        "property": "show",
        "active": true
        },
        {
        "name": "Fight Restriction Zone",
        "property": "show",
        "active": true
        }
    ],
    "display": {
        "category": "airspace",
        "detailedCategory": "Flight Restriction Zone",
        "title": "Manchester/barton",
        "sections": [
        {
            "iconUrl": "https://domain.com//v1/map/icon?icon=warning.png",
            "title": "Flight Restriction Zone",
            "text": "Flight Restriction Zones surround certain protected aerodromes. To Fly in this area you will need permission from Air Traffic Control or the Flight Information Service at this aerodrome."
        },
        {
            "iconUrl": "https://domain.com//v1/map/icon?icon=warning.png",
            "title": "Summary",
            "text": "Red zones are regulated high-risk areas."
        }
        ],
        "actions": []
    }
    },
    "type": "Feature"
},

I wish to loop through each item and apply the style settings to each polygon, based on the returned feature.properties.fillColor and feature.properties.fillOpacity and so on.

Here's my code.

The polygons, markers and lines all display correctly but they are all in the default Blue color and seem to be ignoring my requests to use the fillColor value in the GeoJSON.

// Create an empty layergroup for the data
var LayerAA = L.layerGroup();

// Build Layer: AA
function getAAdata() {
    LayerAA.clearLayers();
    fetch("/v2/mapdata/geojson?n=" + map.getBounds().getNorth() + "&e=" + map.getBounds().getEast() + "&s=" + map.getBounds().getSouth() + "&w=" + map.getBounds().getWest())
      .then(function (responseAA) { return responseAA.json() })
      .then(function (dataAA) {

          // Create a L.GeoJSON out of the data
          var locationsAA = L.geoJson(dataAA, {

              onEachFeature: function (feature, layer) {
                  layer.bindPopup(feature.properties.name, { color: feature.properties.fillColor });
              },

              // Filter out CTA & TMA airspace types
              filter: function (feature, layer) {
                  return (feature.properties.airspaceType !== "CTA" && feature.properties.airspaceType !== "TMA");
              },

          });

          // Add the L.GeoJSON instance to the empty layergroup
          LayerAA.addLayer(locationsAA).addTo(map);

      });

}

Why does this not work?


texts: Capturing GPS coordinates on the ground in Boliva?

texts: Assuming no cellular or internet connection, what is the best way to capture GPS coordinates on the ground in Boliva? Preferably a free approach.  For example, are their any android apps that can do this without a cellular connection?


texts: Finding minimum distance from list of coordinates

texts:
My understanding is that this should work like this, the code will
 start with test list sublist [1,2,3] and move to the points list until
 all lines are run and then move back to test list sublist [2,5,2] and
 continue.  My goal is to obtain a list which will read like this;
min_dist [[1,0],[2,0],[3,0]] because the minimum distance between
 those points are of course the same points.

Am I doing something wrong here?

    import math as m
    test = [[1, 2,3], [2, 5,2],[3,6,8]]
    points = [[1, 2,3],[2,5,2],[3,6,8]]
    distance = []
    min_dist = []
    for i in test:
        x1 = i[1]
        y1 = i[2]
        print i
        for j in points:
            x2 = j[1]
            y2 = j[2]
            ID = j[0]
            print j
            dist = m.sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1))
            distance.append([ID,dist])
            min_dist.append(min(distance))
            del distance[:]


texts: I want to load a raster base layer in leaflet because the base layer usually needs an internet connection I want to load the raster offline .


texts: Serving ESRI layer as WMS through QGIS server

texts:
I added ESRI's World Hillshade layer through Add Layer -> Add ArcGIS Mapserver Layer... and it shows up correctly in QGIS. Now I want to use it as base layer for a online map. I try to offer it as WMS service using QGIS server to later visualize it with OpenLayers.

The problem is that GetCapabilities doesn't include the layer and so I get a LayerNotDefined ServiceException when calling the WMS service.

The only suspicious thing I see on QGIS site is that the layer properties show a exportTilesAllowed=false. Could this be the reason that it's not offered as WMS data, does QGIS server have in mind this property?

Should it generally be possible to serve ArcGIS Map Server layers with QGIS server to web sites?


texts:
I've got a polygon shapefile with 11 features (11 Italian cities) and I want them to merge them into one, however I cannot activate the "Merge selected features" button.

Any suggestion about the reason?

Mind that the project SR and the layer SR are the same and I already verified the topology.

I attach a snapshot of the viewer and a link to download the shapefile.

Download the shapefile


texts: Is there an easy way to snap lines to nearest point of either a line or a polygon? I don't see any snapping options for line layers as suggested on How to merge lines with slightly different endpoints in QGIS?


texts:
I have a Map Package of ArcGIS and my licence expired. I need to use this file now.  I am using QGIS.

How do I do this?


texts: Why do I have various *.sde files on my server?

texts:
I have an old pc at work which has been used to serve a large aerial Photo. I can only imagine that this was done due to lack of space. I have now found 5 other sde files on the server. I have no idea how to read these files, and don´t know what they actually do. I always use direct connect so I imagine this has something to do with a 3 tier connection.

Can anyone tell my what sde files do and why they might possibly be saved on this PC?


texts: How to use the same ROI for multiple LANDSAT images in ENVI?

texts:
Suppose that I have 3 landsat images with the same datum ( WGS-84 ) and projection system ( UTM, Zone 39N ) as follows:

TM LANDSAT 5 image for July 1987
ETM+ LANDSAT 7 image for July 2001
OLI_TIRS LANDSAT 8 image for July 2013

I have drawn a ROI in TM LANDSAT 5 image for July 1987 and saved it in a .roi file and now I want to use that ROI to subset data via roi for all of the images. But I encounter the following problem:

I have opened all of the images in ENVI in multiple displays, restored the roi and then used the option subset data via roi
when I do this for TM LANDSAT 5 image for July 1987, there's no problem and the action is done completely and I can store subset data to a file but for the remaining two image, I encounter the ENVI error:

There are no available ROIs associated with this input file

What can I do?
Because I'm doing a change detection project, I don't want to draw approximately identical ROI's in multiple images???
And also I think regarding the same coordinate system for all of the images, I should be able to draw a ROI in one image and then use it in other images???
Is that true or maybe I'm thinking in a wrong way???


texts:
The Kartograph package looks very promising, so I have been trying play around with it.  As it turns out, I have been able to get the most basic map to work (the "world.json" example), but I could not get any of the shapefiles I actually needed to work. (I am primarily concerned with Washington DC, at least professionally.)  They all throw the same exception:

from kartograph import Kartograph
K = Kartograph()

nbhds={"layers": [{"id":"mylayer",
                   "src": "/home/choct155/ora_public/dc_shp/CompPlanPlanningAreas.shp"
                   }],
       "proj":{"id":"mercator"
               }
       }

K.generate(nbhds, outfile='nbhds.svg')

RuntimeError: projection not named

I tried to follow the world example exactly (changing only the path of the shapefile) and played around with various configurations I saw others use.  Further, as you can see above, I also tried explicitly setting the projection.  Searching this error returns exactly one hit, identifying the LCC projection as the culprit.  All DC shapefiles use the NAD83 datum and Lambert Conformal Conic projection.  Even if I build up from Census shapefiles, as soon as I interact with a DC shapefile (when, for example, mapping blocks to DC geographies) this issue crops up.

I have been able to circumvent this by converting the base file to a different projection, but LCC is really the best projection for my area (and the addition to the workflow is a bit annoying).  Does anyone know of a way around this?

Also, I have generally found less resources than I would expect for some pretty common spatial operations.  In the interest of providing another resource for those learning to deal with spatial data in Python, here is a rough workflow that may be of some use for converting shapefiles from one projection to another.


texts: How to correct border between 2 polygons in R?

texts:
What algorithm to use to correct (average) the borderline between 2 polygons? All I know is that they have their own borderlines very close to each other. Sometimes they may overlap, otherwise not.

As in the image it's visible the common line is fatter than outside borders so it's not ok.
I'd require some average line or just common border - make them adjacent. Those differences are very small so averaging is ok.
EDIT
Given p1 and p2 as SpatialPolygons I tried the following:
sf_p1 = sf::st_as_sf(p1)
sf_p2 = sf::st_as_sf(p2)
sf_u = sf::st_union(sf_p1, sf_p2)
d = gDistance(p1, p2)

The plot remains the same. d takes 3.03 so I guess there's a little gap between them. I'd probably use buffers? But how?


texts: Merging multiple feature classes in file geodatabase with ogr2ogr?

texts:
Is there a way to merge multiple feature classes into one feature class in a geodatabase with ogr2ogr? If this isn't possible with ogr2ogr, is there another method? I'm using the ESRI File Geodatabase driver for GDAL.

If the solution could merge all feature classes in a geodatabase matching a particular string, that would be even better!

update

org2org command I'm using to merge geodatabases, which works well but creates an unwieldy number of feature classes:

ogr2ogr -f 'FileGDB' -update -append merged.gdb gdbtomerge.gdb

Here is the structure of the geodatabase feature classes i'd like to merge:

geodatabase #1 feature classes (multipolygon) to merge:

Sec01_Frm01
Sec01_Frm02
Sec01_Frm03
etc...

geodatabase #2 feature classes (multipolygon) to merge:

Sec02_Frm01
Sec02_Frm02
Sec02_Frm03
etc...

There are no attributes for any of the features except for SHAPE_Length and SHAPE_Area for all feature classes. Ideally I'd like to give each feature a text attribute with the feature class name, but that is out of the scope of this question I think. Just merging all of the feature classes that either match the string 'Frm' or even the same geometry type would work.

update

I think I've gotten a little closer using SQL expressions with ogrinfo:

-drop a table (feature class):

ogrinfo -dialect FileGDB -sql "drop table tabletodrop" geodbname.gdb

-I tried this to merge two tables but I get: '(Cannot acquire a lock.)' with this command:

ogrinfo -dialect FileGDB -sql "INSERT INTO table01 SELECT * FROM table02" geodbname.gdb

Any thoughts on the lock issue or a better way to do this?


texts: Calculating areas in continuous raster of flood map?

texts:
I generated a flood hazard map (raster) of a certain area using ArcGIS 10.5 and classify the result into

low    0.1 ÷ 0.5
medium 0.5 ÷ 1.5
high   > 1.5

After that, I overlay a Land Use map (polygon) in order to identify affected areas. Now, I want to count the area occupied by each class with respect to the land use map. Can anyone suggest methods for that?


texts:
I am very new to QGIS for the primary purpose of creating a slopes layer so we know where slopes are above 15% or 20% or whatever we finally decide.

I have used our DEM to create a slope geotif and have used that to try and create a classified layer, but that is where my problems seem to begin. I used this formula to create my classified layer ((slope@1 < 15) OR (slope@1 > 60)) = 0 and while I get a new layer I am not able to view it, no matter how I manipulate the properites. I also tried to follow the instructions for the display that says to do this:

  9.4.16. Setting the style for the reclassified layers

  Open the Style tab in the layer’s Properties dialog as usual.
  Under the heading Load min / max values from band, select the Actual    (slower) radio button.
  Click the Load button.

  The Custom min / max values fields should now populate with 0 and 1,
  respectively. (If they do not, then there was a mistake with your
  reclassification of the data, and you will need to go over that part
  again.)

  Under the heading Contrast enhancement, set the Current dropdown list    to Stretch To MinMax.
  Click OK.

I do not come up with the 0 and 1 that the instructions speak of (I am not using their tutorial data) but it seems that I should be able to see the layer.


texts: Converting lat long coordinates to shapefile x y values

texts:
I'm trying to figure out where a given latitude, longitude coordinate should appear on a particular shape file.  The shape file was generated by the US government and I doubt it has any peculiar definitions (details below).
As an example, a point in Miami, Florida has real lat,long of 25.77,-80.21.  Looking at the contents of the shape file this region is roughly at a y,x coordinate of 292000, 8930000.  I have tried a python based utm conversion of my lat,longs which goes something like this:
 utm.from_latlon(25.77,-80.21)

that returns the following:
 (579216.2480487269, 2850450.638436528, 17, 'R')

The utm converted latitude number looks roughly close but the longitude number is way off.  Is there some way to convert the longitude to an absolute number by adjusting it with the utm zone value of 17 or am I headed down the wrong path?
The shapefile in question is here: https://catalog.data.gov/dataset/usda-rural-development-property-eligibility-sfh-mfh
which has projection file description of:
PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",
  GEOGCS["GCS_WGS_1984",
    DATUM["D_WGS_1984",
      SPHEROID["WGS_1984",6378137.0,298.257223563]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]],
  PROJECTION["Mercator_Auxiliary_Sphere"],
  PARAMETER["False_Easting",0.0],
  PARAMETER["False_Northing",0.0],
  PARAMETER["Central_Meridian",0.0],
  PARAMETER["Standard_Parallel_1",0.0],
  PARAMETER["Auxiliary_Sphere_Type",0.0],
  UNIT["Meter",1.0]]


texts: Inserting new points to desired distance from the existing ones?

texts:
I am using QGIS to map the existing buildings on the map, we want to build more houses but on a distance not less than 5 km from the existing ones, it means that a place where there is no building in 5 km at least should be marked on the map so that we can know where to put the new infrastructures.

Is there any feature/way that can help me to do it easily in QGIS?


texts: Edit XML from geodatabase using Python element tree

texts:
I am trying to disallow Null values for specific fields within a gdb.  Changing field properties using arcpy.AlterField_managment only changes the properties in the Python object and not in the gdb schema.  My workaround is to export an XML workspace document, make the changes using the element tree module in Python, then import the XML back into the gdb.  When I go to import the XML back into the gdb ArcGIS gives me this error "Invalid XML file".  I removed my code from the equation by just reading the XML workspace document into Python then directly writing it back out.  I get the same error, so the issue is not with my code I am using to edit the input XML.  Also I tried ET.parse() on the output_xml_file and it parses just fine.  So the issue must be that element tree is changing my Esri workspace document when it parses it so that when I try to import back into gdb it doesn't recognize it.

import xml.etree.ElementTree as ET
ET.register_namespace('esri', "http://www.esri.com/schemas/ArcGIS/10.6")
ET.register_namespace('xsi', "http://www.w3.org/2001/XMLSchema-instance")
ET.register_namespace('xs', "http://www.w3.org/2001/XMLSchema")
tree = ET.parse(esri_workspace_doc_xml_file)
tree.write(output_xml_file)


texts:
I'm trying to generate equidistant points on a map, I used QGis to do it, but the problem is that some points are out of the map area. I'd like to consider only the points inside the map.

This is an image with the issue described above:


texts: Synching Two ArcMap JavaScript webmaps to respond to panning

texts:
I am looking for some advice on how to approach a problem I am having. I am looking to build an ArcGIS webmap using JavaScript and HTML 5 that reacts to a user panning on two separate, synchronous maps.

Modest Maps provides an excellent example of Loading Two Maps and Synching Their Panning that helps visualize what I need to do. However, I need to accomplish the Loading Two Maps and Synching Their Panning with ArcMap products.

The Overview dijit in Arcmap seems to be capturing the essence of what I am trying accomplishing. However, the overview dijit  in the widget overviewmap seems to react rather then pan synchronously.

The dijit css feeds into this map and can be used as a resource if needed. I couldn't find anything in it that would accomplish synchronous panning.

I am open to any ArcGIS suite methodology to accomplish the same as the Modest Maps look and feel.


texts: ST_Line_Locate_Point between two spatial relations with PostGIS

texts:
I am trying to calculate the closest edge and its offset for a relation of points.
If using only one point that approach works fine, but if I have several points in the relation I do not know how to do that, because it always return all the edges within the range i specify (here 100). By setting the limit this does not give me the correct answer. What i need is the following:
 join each tuple of the vertex relation with each tuple of edge relation and return the closest edge to the current vertex.

Any idea how to solve this problem?

thanks for any suggestion.

Markus

SELECT L.ID, ST_Line_Locate_Point(L.EDGE_GEO,L.POINT_GEO)*L.LENGTH AS OFFSET,   l.node_id
FROM (
SELECT E.ID, E.GEOMETRY EDGE_GEO, E.LENGTH, N.GEOMETRY POINT_GEO, n.id node_id
FROM
BZ_EDGES E, BZ_VERTEX N
WHERE E.EDGE_MODE=0
AND st_dwithin(E.GEOMETRY,N.GEOMETRY, 100)
ORDER BY node_id, st_distance(E.GEOMETRY,N.GEOMETRY)) L

In Oracle Spatial this approach works fine:

SELECT id, dist, offset, node_id
FROM (
SELECT sdo_nn_distance(1) dist,
SDO_LRS.get_measure(SDO_LRS.PROJECT_PT(
SDO_LRS.CONVERT_TO_LRS_GEOM(e.geometry), v.geometry)) offset,
e.*, v.id node_id
FROM BZ_EDGES E,BZ_VERTEX V
WHERE
sdo_nn(v.geometry,e.geometry, 'SDO_NUM_RES=1', 1 )='TRUE')


texts: Focal statistics using geodesic distance in ArcGIS Desktop?

texts:
I need to compute focal statistics (sum, mean etc.) for raster data spanning the whole globe. I want to use a constant radius e.g. 100km around each cell.

However the data are geographic (WGS84), so setting a radius in map units will result in the geodesic size of this radius changing with latitude.

I've written a little python script to crunch the numbers by re-calculating the filter width with every latitude to keep the geodesic radius constant, but this is slow & a bit inelegant.

Can anyone suggest a better way of doing this in ArcGIS Desktop?


texts: Dynamically load and display shapefiles from PostGIS 2.0 into Mapguide using .NET

texts:
I have a C# mechanism for Loading shapefiles into PostreSQL/PostGIS now I would like to dynamically load it on Mapguide fusion template.In either C# or VB.NET code behind.

Use case is:
1.use uploads a zip folder with shapefile
2. .NET code behind checks for certain attributes then uploads it to PostgreSQL/PostGIS Database
3. upon successful upload, display the layer into Mapguide Fusion template
4.Overlay it over base layers


texts:
I'm searching for an ArcMap (10.2) Tool that merges multiple polygons over a given distance (in the end term I want to select all polygons with an area greater than 7ha). To illustrate my use case here is a screenshot of my problem:

As you can see I selected three polygons that the tool should merge into one polygon. To be more precise I want to merge every feature that has a distance less than or equal to sqrt(100^2+100^2) (this is because the size of a a square is 100x100 meters).
What I tried so far is the Aggregate Polygons tool within the Cartography Toolbox, but the result was unsatisfying:

The red circles show the side effects of filling some gaps I rather want to be ontouched. Also the area indicated by the green circle is not a part of the merged polygon.
The following is the generated python snippet of my tool execution:

arcpy.AggregatePolygons_cartography(in_features=...,out_feature_class=...,aggregation_distance="142 Meters",minimum_area="0 SquareMeters",minimum_hole_size="0 SquareMeters",orthogonality_option="ORTHOGONAL",barrier_features="#")

Am I using the wrong tool for the job?

My current workaround is to copy the featureclass; create a buffer within the copy; dissolve (with multipart none); select the features with area > 7ha; select by intersect; copy the selection to a new feature class

After talking this over with a tech from our group I'm afraid the data is not given as a raster file, but only as a feature class (and converting it back and forth again wouldn't be ideal)

For clarification, here is a simplified problem of three diagonaly adjacent 100x100 meter polygons(/features @I). @II/III, we're searching for a solution that scans for neighbors not only in a horizontal and vertical way, but also in diagonal directions. e.g. @IV scans for neighborcandidates - @V finds one, adds the square to the feature(@VI) and scans for new neighbors. The result should look like @VII with all three squares merged together and the according data records should state a shape area of 300 sq m.


texts: I have successfully published a Geoprocessing Tool as GP services, I had called it with success, but when I try to call it, now it gives me this error GPServer does not exist or is inaccessible for no reason. MapServer and GeocodeServer work fine I am having a problem just with the GPServer.


texts:
It seems I am having a problem with ST_Split function. I am trying the following but it does not work. Can anyone spot any problems with the syntax?

select ST_Split(
   ST_GeomFromText('LINESTRING(-3.164154576738156 55.926335330437809,-3.16323392764116 55.926283382919387)'),
   ST_GeomFromText('LINESTRING(-3.16348554978961 55.92654456162947,-3.163507914241407 55.926112217812459)'));

I have tried some more SQL (including the examples in the PostGIS documentation) with st_split and it does not seem to work...


texts:
when i tried to georeference a printout from map.geo.admin.ch, either the result got turned some degrees, or the height got compressed or both.

using the helmert transformation gave me the best results. the result was not compressed - only rotated.

it then looks like:

the coordinates for the reference points i userd were taken from the coordinate grid on the website. i used the four corners.
top left: 46.86851  9.99669
bottom left: 46.76063   9.99157
bottom right: 46.75528  10.22048
top right: 46.86315 10.22609

i tried it with a screenshot of google maps with the same results...
is there a way to georeference that print without rotating it?


texts:
I need to be able to walk through many mxds to replace text strings but able to skip a character place. For example, replace:

'area %d'

with

'area %d - this is some other text'

output:

'area 2'

'area 2 - this is some example text'

'area 3'

'area 3 - this is some example text'

This is what I have so far. It works for replace a specific string but I have not been able to figure out the placeholder variable.

import os
import arcpy

mxd_path = r'C:\Users\username\Desktop\test'

for dirpath, dirnames, filenames in os.walk(mxd_path):
    for filename in filenames:
        fullPath = os.path.join(dirpath, filename)
        basename, extension = os.path.splitext(filename)
        if extension.lower() == ".mxd":
            has_data_source = False
            mxd = arcpy.mapping.MapDocument(fullPath)
            for elm in arcpy.mapping.ListLayoutElements(mxd, 'TEXT_ELEMENT'):
                elm.text = elm.text.replace('area %d',
                                            'area %d - this is some other text')

            mxd.save()
            del mxd


texts: Looking for GIS server for Pictometry that can do Bing-like "birds eye" tilted view and rotation

texts:
I need an open source GIS server that can present these Pictometry images, I guess they are also called 2.5-D. The camera on a plane takes 5 shots, one orthogonal, other 4 in all directions. North, west...

2.5-D images are the main focus here, user would have to be able to rotate the tilted view in 4 directions. Additionally I would like to have those standard GIS server features like vector layers, queries, labels, redlining...

So far I have decent experience with MapGuide, open source and the the Autodesk one, and I can't recall those features. No rotation or tilting at all, I might be wrong. I have very basic XP with GeoServer and MapServer.

I'm looking for advice in which direction to go with research. I was thinking GeoServer is the best bet, but what interface...

Also, I can't really say that I've been working with any kind of 3-D desktop GIS. I would take some advice on that too.

Any suggestions?