texts: Digitizing contour lines from a picture in QGIS

texts: I have a map with isobaths (water depth), like in this picture, and I would like to digitize them. Is there any tool or procedure in QGIS that can recognize the lines in the picture and create a vector file automatically?


texts: Normalization where the smallest grid value gets the biggest value with “rescale by function”

texts:
I tried to normalize a density analysis between 0 and 100 with the tool “rescale by function”. The tricky thing is now, that the smallest value out of the density analysis should get the value 100 (here: 0) and the biggest value should get the value 0 (here: 5.28). I used the tool “rescale by function”, as shown below.

But I am struggling with the result, because the lowest grid value is negative, but it should be 0 (null). Does anyone have a solution for this problem? The result grid should have the values 0 – 100, where 100 was the smallest value from the density analysis and 0 the biggest value from the density output grid.
I am working with ArcGIS 10.6.


texts: Differences between 'target' and 'selected' properties of an event in OpenLayers

texts:
In the case of a browser event in OpenLayers, what exactly is the difference between:
select.on('select', function (event) {
    let selectedFeature = event.selected[0];
}

and:
select.on('select', function (event) {
    let selectedFeature = event.target.getFeatures().item(0);
}

select being an ol.interaction.Select instance.
Here's an output of such event on the console:

They seem to give the same results at first glance.
Underlying this question, I also naturally wonder if one is better than the other, or not recommended.

As a background context, I first tried for hours to figure out why the following didn't work before finding that it did work with target instead of selected:
select.on('select', function (event) {
    let selectedFeature = event.selected.getFeatures().item(0);
}

From a natural language point of view, this latter makes the most sense to me.
I have found some info on the target (but I do not know if this relates to the same thing than in OpenLayers) at this URL; https://developer.mozilla.org/en-US/docs/Web/API/Event/target but I cannot find anything on the selected method/property.


texts:
we serve our Geographic data via our flagship application Dorset Explorer. This uses OpenLayers and GeoServer to serve our Geographic data, including OSM and OS Basemaps and a mass of other data such as Aerial Photography, Schools etc.

My question is, what sort of Server Specs should we be using for GeoServer. We will only have available to us one server for delivery of tiles via WMS (and perhaps in the future we may be using some WFS on it as well). We get between 100 and 200 users per day, but this will hopefully start rising as soon as we are confident that the servers will hold up! Is it just a matter of stick in the highest specs we can afford or are there some limits or theoretical maximums.

Any pointers would be much appreciated!


texts: Geopandas - Import CSV: make polygons from X,Y Coordinates based on ID

texts:
I'm working on a script, where I want to import a CSV file with coordinates and an ID which should be used to make a polygon shapefile from those.

CSV file looks like this (example):

ID, X, Y

10,15.116686,61.483157

10,17.114749,62.483098

10,17.113456,62.492142

11,14.123456,61.123456

12,12.345678,61.123456
.....
Etc.

So far I managed to make a point output using only the coordinates, but not the ID:

import pandas as pd
from geopandas import GeoDataFrame as gdf
from shapely.geometry import Point, Polygon

temp = r'C:\Temp\%s.csv'
df = pd.read_csv(temp % 'csvfile')
geometry = [Point(xy) for xy in zip(df.x, df.y)]
geo_df.to_file(filename = 'datatest2.shp', driver ='ESRI Shapefile')

This turns out great. But what should i do if i want all the ID's with 10, gather those coordinate points, and then make a polygon out of those coordinates? - and ofc. with ID 11,12,13 etc. as well?


texts: Why do vectors disappears when zooming in?

texts:
I have imported a .gpx file into qgis and separated the track waypoints using Vector > GPS > GPS tools.  When viewing at full extent I can see the original track and component waypoints.  However when I zoom in, the track lines disappear but the waypoints are still there. I have tried:

checking track layer scale dependent rendering is OFF
exporting track layer as .shp and adding to map
moving track layer to top of layer stack
toggling map rendering

but still cannot view track when zooming in.  Any ideas?
Thanks,

Tony

[EDIT] I have uploaded Pic1 which shows zoomed to the layer extent ('GPX tracklog export').  Pic2 is when it has zoomed in with mouse wheel and part of the shapefile has diappeared. Note these are the gpx file tracklog exported as a shapefile.


texts: Transfer domain descriptions from coded values to string text at a new field

texts: Further to my earlier question,I'm trying to develop a script that eventually should contain three parameters, feature layer, domain field and a new text field. How can i write a loop that will update domain descriptions to a new text field ? I understand i need to use arcpy.da.ListDomains but i have some trouble writing the update loop.Thanks


texts: Can arcpy let the user enter a query expression for select layer by attribute?

texts:
I have written something like this:

inFc = sys.argv[1]
FLayer = sys.argv[2]
FValue = sys.argv[3]
arcpy.MakeFeatureLayer_management(inFc, FLayer) #feature layer created
arcpy.AddMessage ("Layer has been created")
arcpy.SelectLayerByAttribute_management(FLayer, "NEW_SELECTION", " FValue ")

In my parameters, I have inFc set as an input feature class, FLayer as the output feature layer, and FValue as an input query expression. I also clicked "obtain the data from the input feature class" for the FValue query. When I do this, however, it presents an error. I'm sure it has something to do with the user-entered query.


texts:
I have a Date column in my feature class and I want the date to appear in "dd MMM yyyy" format as a label to a polygon. I am using Arc Map 10.1.
I tried taking Jscript.

var Dt  = new Date();
Dt = Dt.toString("dd MMM yyyy");

But it gives error. I can implement the solution in jscript ,vbscript and python if available.


texts: How to extract data for a specific area from OSM?

texts:
I need to extract all streets available on a particular area of a City. I'm able to extract all streets name using this query:

select distinct P2.*
from
planet_osm_polygon P1,
planet_osm_line P2
where P2.name is not null AND P1.osm_id = -80277 and ST_Intersects(P2.way, P1.way)
order by P2.name

(in this casa -324211 is the city of Toronto). Now, I need to extract the streets name of only a particular subarea of the City. I can found the information about the subarea on planet_osm_point table but, how can I extract only the line or rows are in one particular subarea?

With this I can extract all subarea of Toronto:

select  P2.osm_id, P2.boundary, P2.name
from
planet_osm_polygon P1,
planet_osm_point P2
where P2.name is not null AND P1.osm_id = -324211 and ST_Intersects(P2.way, P1.way)  and P2.place = 'suburb'
order by P2.name

is there a way to do that?

--- UPDATE ---

One solution is, starting from a point, create a square or round polygon and extracts all streets that intersect this new polygon, something like that:

select asText(way) from planet_osm_point where osm_id = 349354234

and, when I use the results here:

select  P2.*
from
planet_osm_line P2
where P2.name is not null AND ST_Intersects(P2.way, (ST_Buffer('POINT(-8840304.91218709 5429889.22642508)', 50) ) )

I get this error message: Operation on mixed SRID geometries

On a blog I read how to update my SRID table geometry column but, is it the right way to fix that?

thanks for your help! Any suggestion (not only a solution) will be accepted and very helpful!

Andrea

PS: I imported data using osm2pgsql but, if I understand well (I'm newbie here) I can found more information using osmosis, correct?


texts: Python script for if/elif condition in field calculator

texts:
I need some help in calculating the last column (OD_pairs.Ai_) based on the condition;
if travel_cost has a value below 10, OD-pairs.Ai_ takes the value of the first column(Aij_OjFCij..)
if travel_cost has a value above 10, then return 0 in OD_pairs.Ai.
So basically the third field to be calculated is either the corresponding value in the first field or zero, depending on the value in second field.


texts: Truncate arbitrary File Geodatabase tables in FME

texts:
I have one FGDB with over 30 tables. Also have few WB and they write data to output tables in FGDB, but each of WB write data to different tables but still at the same FGDB. Output tables are being truncated first (set option Parameters). But I also need the rest of tables in this database to be truncated every time WB is started, no matter which one.
Because the problem is that FGDB is going to be updated only with data accoridng to spesific WB, but there are still data in other tables not used by WB.
I hope someone understand my issue.
Tnx.


texts: Calculate stream locations and Strahler order from DEM in QGIS

texts: I am trying to locate streams and determine their order in QGIS. I have tried using several of the procedures described in answers to similar questions, but none of them seem to work. I have filled the sinks in my DEM, and then none of processes produce valid results. They either do not create the layers at all, or when the layers are created there are areas of internal flow that do not exist in the real world. I have GRASS included in the QGIS distribution on my Linux OS, and SAGA on my Windows OS, but neither of them contain both. I would prefer to do this process in either one or the other, but I can utilize tools from both if necessary.


texts: Query not working on Virtual Layer but working on the real layer

texts:
I'm using QGIS 3.10.10 because of project compatibility. I have the following issue with a query used to filter a virtual layer:
If I try to apply a filter to the layer INFRASTRUTTURE as follows:
SELECT SUM(L_TOT) FROM (SELECT ((N_TUBI12 - XNTUBI12)*LUNGHEZZA) AS L_TOT
    FROM INFRASTRUTTURE
    WHERE N_TUBI12 <> '0' or N_TUBI12 <> null)

I get the error:

An error occurred when executing the query. The data provider said:
OGR[3] error 1: SQL Expression Parsing Error: syntax error, unexpected
'(', expecting string or identifier. Occurred around : SELECT
SUM(L_TOT) FROM (SELECT TIPOLOGIA, N_TUBI12, XNTUBI12, L*

If I try to filter in two steps (creating a new layer L_TUBI_12 in the process to store the result), no error is found and it works perfectly:
SELECT ((N_TUBI12 - XNTUBI12)*LUNGHEZZA) AS L_TOT
    FROM INFRASTRUTTURE
    WHERE N_TUBI12 <> '0' or N_TUBI12 <> null

SELECT SUM(L_TOT) FROM L_TUBI_12

Now, if I create a virtual layer and I try the first query, no error is shown but the result is NULL. I tried forcing the output as int (because it was Qstring by default) but the result is still NULL.
SELECT SUM(L_TOT) FROM (SELECT ((N_TUBI12 - XNTUBI12)*LUNGHEZZA) AS L_TOT
    FROM INFRASTRUTTURE
    WHERE N_TUBI12 <> '0' or N_TUBI12 <> null)

The problem seems to be the additional operation - XNTUBI12 because, if removed, it gives a number as result. The subtract operator is basically not working in the virtual layer.
Do you know why this behaviour? Is there any workaround?


texts: English-taught bachelor programs in cartography and GIS science?

texts:
Is there any English-taught bachelor program in Cartography and GIS science in Belgium, the Netherlands, Austria, Germany, Denmark, Sweden, Finland?

I've done massive homework about that, but haven't succeeded. If you know any course or program like that, please let me know.

As the only international language that I know is English, therefore I can only enroll on English-taught programs. Moreover, programs in USA, UK and Canada are unaffordable to me, so I want to do my study in the above countries.


texts: Creating snapshot in ArcMap of layer to use for legend?

texts:
I am trying to create a snapshot using FastStone Capture version 5.3 to be use for the legend because it is diffcult to match their color especially is you are using raster or transparent ?

How do I save the snapshot to file and then use that file for the legend ? Any suggest ?

I am using 9.3.1 SP3.


texts: Cannot connect to current version of enterprise database (PostgreSQL) through PHP

texts:
Very new to PHP, and working on a simple edit page. I'm connecting to multiversioned Postgresql database that I set up in ArcMap using a PHP/HTML file. I'm using the pg_connect command with host/dbname/user/password as the input.

$MyPG = pg_pconnect('host=sv01213 dbname=surveys user=sde password=test', PGSQL_CONNECT_FORCE_NEW);

It seems like when I make edits to the attributes through the PHP interface, I'm not seeing any changes have saved when I view the attributes in ArcMap. And also visa versa, when I make changes to the database through ArcMap, I'm not seeing the changes in a PHP query table. Even when I changed between the two versions of my postgresql database, I'm not seeing the PHP changes.

Does anyone have this same behavior? It almost seems like the PHP connection is making a separate static connection to the database or connecting to a alternative version. Is there a way to specify which version of the database that it finds?


texts: How to specify the lat/long and bounding box of MapServer WMS layer in OpenLayers

texts:
I am using MapServer to add a WMS layer on a basemap using OpenLayers 3. The WMS generates a png image from a shapefile, which I need to overlay on a basemap. Below is the code which I used to overlay the WMS layer.

var zone_wms = new OpenLayers.Layer.WMS("FEMA Flood Zones",
"http://localhost:81/cgi-bin/mapserv.exe?map=/MapServer/maps/femawms.map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=FEMA_NFHL_POLY_STATUS&STYLES=&SRS=EPSG:4326&BBOX=-69730.2,777408,331075,959674&FORMAT=image/png", {
  layers: "Flood_Hazard_Zones",
  transparent: "true",
  format: "image/png"
});

map.addLayer(zone_wms);

However, when I do that, the layer is repeated on each of the tile on the basemap. For example, when I am at a zoom level=1, there are 8 images (one of each tile) displayed. When I am at zoom level=2, there are 32 images displayed on the basemap... and so on. Also, I am not sure how to specify the exact location (lat+long+bounding box) of the WMS layer that is overlaid on the map.

I am sure, you would have figured it out that I am newbie to the GIS world..so please forgive me for my ignorance.


texts: How to Export Shapefile Attribute to Text Layer in CAD Using FME

texts:
How can I export a Shapefile's Attribute to Text Layer in CAD Using FME.

For example I have a Shapefile of Parcels and it contains a ParcelID field. Now I would like to display that field attributes as CAD Text Annotaion layer.
Thanks


texts: OpenLayers WMS - TypeError: this.resolutions is undefined

texts:
I'm trying to setup a map with WMS and WFS Layers by using Geoserver.
However, when I'm combining WFS and WMS (with some Google basemaps) together in one map I got following error in Firebug:
TypeError: this.resolutions is undefined
which is located in the http://openlayers.org/dev/OpenLayers.js JS-file
Then I can see the WFS layer but not the WMS layer, furthermore also navigating in the map becomes very difficult...

The WFS Layer look like this:

       var tracks = new OpenLayers.Layer.Vector("Tracks", {
    strategies: [new OpenLayers.Strategy.BBOX()],
    styleMap: style,
    protocol: new OpenLayers.Protocol.WFS({
        version: "1.1.0",
        url: "http://localhost:8080/geoserver/wfs",
        featureNS: "http://shapes.net/shapes",
        featureType: "tracks",
        geometryName: "the_geom",
    })
});

Well, thise works quite good but as soon as I'm adding also a WMS Layer into the map i got the Error

var corine =  new OpenLayers.Layer.WMS( "OpenLayers WMS",
                "http://vmap0.tiles.osgeo.org/wms/vmap0",
                {layers: 'basic'} );

It doesn't matter if I take a WMS from my Geoserver or from an other external source (like in the code here), I always got the error when I'm activating this layer.

Btw. if the layers are not combined the works fine...

My assumption is that there could be a problem with the projection. Both layers are defined in Geoserver with EPSG:4326.

In the init.js I'm using this projects here:

projection: new OpenLayers.Projection("EPSG:4326"),
displayProjection: new OpenLayers.Projection("EPSG:900913"),

Any ideas?


texts:
I have this web resource XML description. From what I saw, it's ArcGIS API format, but I can't find anything else.

What standard is it, and is it possible to connect to this resource in QGIS?

http://www.metadados.geo.ibge.gov.br/geonetwork_ibge/srv/por/xml.metadata.get?uuid=b3c58b5e-10d5-4026-a4f5-dab1393054a8


texts: Show representation properties in object inspector in ArcMap 10

texts:
When you inspected a feature in the object inspector in ArcMap 9.x (opened by clicking the  "Attributes" button in the editor toolbar), you were able to see the properties of ALL the features representations, as individual tabs in the dialog (behind the attribute tab).

In ArcGIS 10.0, I have not been able to do this. Only, if the representation is active, can I see its properties in the inspector. This has ENORMOUS (!!!) costs to my work.

Until now, my workaround has been to add the layer several times with the different representations, but I am kind of hoping, that there is some way to make. However doing this on a country-wide road network, makes it painstakingly slow to work in ArcMap, so i really hope that one of you has some way of making my problem go away.


texts: Add custom variables to leaflet polyline

texts:
I'm drawing a route on a map with leaflet. I'd like to click on any point in the line and display the speed and time at that point in a popup ($speed[$i]).

latlngs = [<?php for ($i=0;$i<count($latts);$i++){echo '['. $latts[$i] .', '. $longs[$i]     .']';if($i!=count($latts)-1){echo ',';}} ?>],
itinerary = L.polyline(latlngs, {color: 'red'});
itinerary.addEventListener('click dblclick', function(e) {
    console.log(e);
    var res = e.latlng; // that's nice but I'd like to be able to get e.speed, etc...
});
map.addLayer(itinerary);

EDIT: I'm building the polyline with an array of Latitudes and Longitudes. I've got a Speed array of the same length. When the user clicks on a point on the line, I'd like to get the closest recorded Latitude and Longitude, and display the corresponding speed for that point (in a popup or in an other div).
Eventually the user can slide his mouse along the line and see the speed .


texts:
I'm trying to extract administrative relations from a country extract of Ukraine (downloaded from http://download.gisgraphy.com/openstreetmap/pbf/UA.tar.bz2). Ultimately I'm wanting to get the OSM administrative relations converted to GeoJSON.

I'm using Osmosis to parse the PBF and write to a new PBF with just the administrative relations. When I run Osmosis as below, I seem to losing some shapes. In the case of Ukraine I'm losing the relation for the country (OSM ID 60199):

$ osmosis --read-pbf file=UA.pbf --tf accept-relations "boundary=administrative" --used-way --used-node --write-pbf UA_administrative.pbf

Any idea why this leave out some shapes? Oddly enough, if I use Osmosis to insert directly into a PostGIS database with --write-pgsql, the relation for Ukraine is there.


texts: Leaflet: improve code for show/hide legend icon

texts:
I'd like to know how I can improve the code below, making it better and easier to read. The idea is that when the user turns on / off the layer, their representation in the legend is also shown / hidden:

The complete example is here: https://codepen.io/marcellobenigno/pen/MZRZgg

in this part of the code is configured the functionality

L.Control.Layers.include({
        getOverlays: function () {
            var control, layers;
            layers = {};
            control = this;
            control._layers.forEach(function (obj) {
                var layerName;
                if (obj.overlay) {
                    layerName = obj.name;
                    return layers[layerName] = control._map.hasLayer(obj.layer);
                }
            });
            return layers;
        }
    });

var legend = L.control({position: 'bottomright'});

legend.onAdd = function (map) {
    var div = L.DomUtil.create('div', 'info legend');
    var title = '
LEGEND:
';
    div.innerHTML = title;
    div.innerHTML += ' Free Bus';
    div.innerHTML += ' Campus';
    div.innerHTML += ' Bicycle Rental';
    return div;
};
legend.addTo(map);

function showLegend(layer) {
    var item = ".legend > span:contains(" + layer + ")";
    $(item).show();
}

function hideLegend(layer) {
    var item = ".legend > span:contains(" + layer + ")";
    $(item).hide();
}

map.on('overlayadd', function (eventLayer) {
    showLegend(eventLayer.name);
});

map.on('overlayremove', function (eventLayer) {
    hideLegend(eventLayer.name);
});

$.each(control.getOverlays(), function (index, value) {
    if (!value) {
        hideLegend(index);
    }
});


texts:
I want to calculate the area of each polygon in a shapefile.

How do I control what the units of measurement are?

def getArea(shapefile):
    """Print area of polygons"""

    Rows = ARCPY.SearchCursor( shapefile )

    for row in Rows:

        tmp_area = row.area # this assumes a field name 'area' exists
        print str(tmp_area)

I looked to see if there is any setting in the arcpy.env. module but could not find one. If I do set the units in my python script, does that mean that all aspects of the shapefile remain unchanged? For example, if the shapefile had units of feet, in my script I change the units to KM, will the shapefile still have units of feet afterwards?


texts:
I have a Shapefile imported into a QGIS project as layer. And I have created a Style for that layer using a Column from the attribute table. To name it I have grouped several landuse-classes to a group (for example: crop, corn etc. would be grouped to arable-field) to give them a single color in the main window. The names for the landuse classes are in the column named "landUse" in the table.

What I want to do is export to a Raster file (ESRI ASCII Raster including a raster attribute table preferably) but using my grouped classes instead of the normal landUse field in the attribute table as cell info.

When I click on the "vector to raster" menu item in the "raster" menu I can only choose a column from the attribute table.

One of my ideas would be creating a column in the attribute table and giving it a value depending on the value in the column "landUse" and then use this column to rasterize. How do I do that and would that even work or are there other ways?

I am not experienced in GIS. I've tried if a few examples/tutorials similar to my problem to see if they could be applied, but yet nothing works.


texts: How to get the a FeatureLayer full extent?

texts: Is there a way to get the extent of a FeatureLayer (with all its features, not only the ones on the map)?


texts: Obtain Polygon Extent with Multiple Polygons within Shapefile

texts:
I have a shapefile that contains multiple polygons. I would like to extract the extent of each of those polygons within that shapefile. The majority of answers that I have been able to find use arcpy and I do not have AcrMap or ArcPro, so I cannot use it. So far I have code that will obtain the extent of the shapefile, but not each individual polygon within the shapefile. I am using fiona and shapely, but am open to use gdal or other FREE python libraries.
from shapely.geometry import shape
import fiona

path = 'D:/data/shapefile/'
file = '__932.shp'

c = fiona.open(path + file)
field = next(iter(c)
extent = shape(field['geometry']).bounds
print(extent)

# (465452.3213890702, 4924992.33678632, 465453.50800855865, 4924993.2240787465)


texts: how to join a road network into as few lines as possible?

texts:
I am processing a polyline file that contains multiple lines representing city roads. I would like to join the lines together into as few lines as possible. I've tried using the 'join multiple lines' plugin, but that creates additional lines that I do not want. Is there a way to only join all the lines that have the exact same endpoint/startpoint?

For context, I'm printing the vector file to a laser engraver. Unfortunately the laser engraver draws the lines in a (seemingly) random order, resulting in a very long print time (hours). I am trying to join the lines together to force the printer to draw connected lines in a logical order. I've tried the 'merge' command, but this hasn't worked because it doesn't change the order of the actual nodes.

I am using QGIS but am open to using other open source software to solve this problem.

Figure 1: Original roads in individual segments

Figure 2: Roads after using 'Join Multiple Lines' plugin


texts: QGIS: Create grid - stuck on 'Processing algorithm...'

texts:
Trying to create a grid from a shapefile (using canvas extent), but when I click 'run', it tells me it is 'Processing algorithm...' and then hangs on this indefinitely without progressing (with my Mac often telling me to Force Quit QGIS as the programme is not responding). EPSG is 4277 - OSGB 1936.

I've tried this on QGIS 3.4.1 and 2.18.0 and the same is happening on both.

Am I doing something wrong?

I'm new to using QGIS (following someone else's instruction sheet) and this is the extent of the instructions I have been given.

What should I be selecting for the different options in the 'create grid' interface?

This is what I've got so far. I know those warning triangles shouldn't be there, but I haven't been able to work out how to get rid of them. Increasing the degrees figure doesn't seem to change anything.


texts:
How can I temporary disable zoming/draging the Mapview in Leaflet.js
Tried so many ways but without any luck.
It's important to make it temporary and I also need the option to enable again.


texts: Imagery Alpha and Mask channels: when should each be used?

texts:
Gdal supports both Mask and Alpha channels. In what circumstances are they best applied? How to choose between them?
For example the nearblack utility to change image collar pixels to nodata can use both -setalpha and -setmask for this purpose. They seem to be functionally identical. Does it matter which is used? Also both can be used at the same time. Does this accidentally arrive from how it's implemented or is there a purpose for using simultaneously?


texts: How to calculate the area between two contour line?

texts: I created contours by using SRTM30 data. Now I need to find area between two contours like in between 0 m to 50 m depth and 50 m to 100 m.


texts: Rasterize point layer with weighted distance

texts:
I have a point layer representing wells in a region. i am able to rasterize it and write the well density into the raster as integer, so far so good.

What I aim for is to create a raster with values in each cell that weight the distance to well-bearing-cells. Raster cells with wells in them (cell value > 1) should get assigned a value of 1.0, cells more than 100km away from well-cells should get a value of 2.0. The cells in between should get values ranging from 1.01 to 1.99, depending on their distance to a well-bearing-cell.

I am not the most profound GIS profesisonal and lack a good approach to that problem. The heat map function seems to go into that direction but i dont get into the finer controls of it. The "weigthed distance to power" tool could be it but if the regression is to a power than its a exponential function (I might be wrong here?) while I am looking to implement a linear function.

Hope someone might have an idea.

I am used to QGIS but getting into R, so I might be able to handle solutions going in that directions.

schematic abstraction of how the resulting raster should be composed. Cells with wells get value=0:


texts: Best way to exclude 0,0 from geography results (PostGIS, GeoDjango)

texts:
Using PostGis and GeoDjango with a model having a Geography field I'm trying to exclude results where the location is 0,0. I get the error PostGIS geography does not support the "exact" lookup. The only thing I can think to do is exclude points that are 0 meters from 0,0:

.exclude(location__distance_lte=(Point(0, 0), D(m=1)))

Is there a better way?


texts: Adjusting / snapping a polygon to another in ArcGIS Desktop

texts: I have a feature class and/or Shapefile that consists of 2 complex polygons.   I would like to remove the small gaps between these polygons 1 & 2,  but I want  polygon 1 to remain exactly the same and only move polygon 2.   I've tried using topology, but I don't see how to make a rule that keeps a record the same and only moves another.


texts: Using Label classes vs Definition queries - which performs the best?

texts: Which would provide the best performance for working within ArcMap with very large datasets stored in an ArcSDE geodatabase? Multiple label classes using complex SQL queries to label specific features, or multiple copies of the same layer using the same complex SQL queries set as a definition query for each layer?


texts: Add new feature to QGIS layer with pre-defined attributes using pyqgis

texts:
I am working on a pyqgis script where I allow the user to add features using the addfeature map tool, where the user can create the geometry of the feature by clicking on the canvas, and after the geometry has been created, the user will be prompt up with a table where he can enter the attributes of the feature.

The question is, some attributes of the feature should be pre-defined, and some can be entered by the user. I found out that by using following code I can suppress the prompt.

layer.editFormConfig().setSuppress(QgsEditFormConfig.SuppressOff)

But I still want the user to be able to enter most of the attributes after the new feature has been added, so it is not a good solution.

Is there another way where I can define some attributes after the feature  has been added, and then prompt the table? In that way, the pre-defined value should be shown in the table when it is opened, so the user does not have to give value to them again.


texts:
I was wondering if there is a way to have CartoDB's torque read end dates for point files, the same way that ArcGIS desktop can. I am trying to plot the foundation of certain cities, and then remove from the map those that were subsequently abandoned. My data would look something like this (assuming the end date of the visualization was the year 1650)

PLACE        FOUNDED     ABANDONED
City_1       1585        1650
City_2       1600        1625
City_3       1615        1620
City_4       1630        1650

As it stands, I have the choice of (a) having the cities appear/disappear according to their foundation date or (b) having all of the cities accumulate on the map. Neither of these solutions work because for a given year (i.e. 1622), I will either show cities that no longer existed or hide some that continued to be there. I know that there is a similar discussion on this matter here, but the solution of adding a row for each year is not feasible in my case.


texts: Is it possible to make a copy of legend in ArcGIS without losing properties?

texts:
I want to make duplicates of my legend in ArcGIS 10.0 that have identical formatting/properties.  Is there a way to copy+paste while preserving properties or a way to save legend properties as a model to apply (e.g. width between columns = 30; font size for layers=12 pt, etc)? Also, preserving the capacity to update giving map layers present and such functionality that is lost when converting to graphic.

For example, instead of preserving the patch width at 30 pt and the height at 15 pt, the patch size in the copied legend is width=290 pt and height=5 pt. Many other properties are also skewed (spacing between, font size for layer/heading/etc, ...).

Context: I am working with a map in ArcGIS 10.0 for which I will be exporting multiple layouts as jpeg/pdf files.  I would like to make copies of my legend to keep in the margins of the data layout view to keep track of changes for different maps.  But when I right-click legend "copy" and then click elsewhere "paste", I lose ALL the formatting for my legend.


texts: How to display table with more than one image from postgis in mapserver?

texts:
I have 3 tif files in one postgis(2.1.5 ver.) raster table. ( load with raster2pgsql)
And i have mapserver MapServer 6.1-dev.

I need to display table from postgis in qgis via mapserver. P.S I trying to do this like with 1 image in raster table but qgis show me just white screen :( Any suggestions?


texts: only 1 tile visible in GE after gdal_translate -of KMLSUPEROVERLAY

texts:
I have a 3-band ASTER scene in any of several formats (PCIDISK, GeoTIFF, Imagine), already reprojected to the geographic coordinate system. I can view and manipulate the scene successfully in PCI Geomatica and ESRI ArcMap. I want to create a KMZ version of it, and following in part this other GIS StackExchange post, I tried using:

gdal_translate -of KMLSUPEROVERLAY myscene.tif myscene.kmz

but the result is invisible in Google Earth (although it zooms to the correct area). It remains invisible when zooming way out and way in, except that the lower-left tile of the 16 tiles it created does show up when zooming in. It, however, is completely black as it is in the no-data collar around the image data.

Also, if I unzip the KMZ and open any of the individual tile KMLs that were created, they open just fine, and the ones that are not in the collar show the correct image data. So I suspect there's something wrong in the way the output KMZ was formed (which uses Regions and NetworkLinks, KML elements I have not worked with before). But I don't know whether this is due to input data or format problems, to incorrect or missing parameters to gdal_translate, or something else. I've tried several input formats with the same results each time.

Tips and suggestions (if not outright answers!) much appreciated!


texts: Adjust font size in QGIS Open Layers Plugin

texts:
I have a problem with the export of open layers maps via the printcomposer. In the printcomposer window the font size of the map labels is easy readable, but if I export the map as pdf or jpg the label size decreases so that they are not readable any more. Decreasing the export resolution leads to a bigger font size in the export file, but also to a very poor quality of the pdf or jpg.
What I need is a high resolution export with a readable font size.

Does anybody know how to solve this problem?


texts: Retaining custom CSS parts when editing using CartoDB GUI?

texts:
I have set a bunch of zoom-level settings in a cartdb map layer, i change opacity, line thickness, etc.

however, when i go into the web GUI and do something simple like change the colour of a line, it blows away all of the custom CSS code.

Is there any way to continue using the GUI after writing additional CSS lines without having the manually added CSS get stomped on?


texts: Extracting outtermost points from a irregular grid of points

texts:
Good day to all,

I am using Arcmap 9.3, and would like to remove all points inside an irregularly placed grid, except the outermost (boundary) points. For example, I would like to keep the black dots, while removing the gray dots. I have over a million points, therefore manually editing the shapefile won't be very efficient. Can I get any pointers/suggestions?

Regards
Kurt


texts: R visualization of the correlation of a point distribution to the attributes of the polygons they are located in

texts:
Situation:
I have two spatial data frames, one of type polygon (areas with different average wind speed levels) and one of type point (wind mills).
I read the polygon shapefile with readOGR, the points were imported as csv and I added coordinates and a projection to it.
What I have:
I use the function over() from the sp-package in combination with summary to determine how many wind mills are placed in areas of each category of average wind speed. Afterwards, I can compare this distribution to the general distribution of wind speed levels and get some kind of result.
Question:
I want a plot that visualizes the relationship between the positions of the wind mills and the wind speed levels. (So hopefully see that wind mills are preferably placed in polygons with high wind speed.) How can I accomplish this?

library('raster')
library('rgeos')
library('sp')
library('rgdal')

polygons = readOGR(dsn = "./North_South_Dakota_Wind_High_Resolution", layer = "ndsd_50mwind")
wind_mills = read.csv('NDGISHubData/NDHUB_WINDTURBINES.csv')
# Assignment modified according
coordinates(wind_mills) <- ~ LONGITUDE + LATITUDE
# Set the projection of the SpatialPointsDataFrame using the projection of the shapefile
proj4string(wind_mills) <- proj4string(polygons)
# get wind mill distribution
windpower_values = over(wind_mills, polygons)
# remove NAs
windpower_distribution = lapply(windpower_values, function(x) x[!is.na(x)])
# create summaries for correlation comparison
# WPC is the attribute containing the wind speed level
a = summary(windpower_values$WPC)
b = summary(polygons@data$WPC)
# ToDo: spatial analysis?


texts: Creating a hyperlink within a popup using Leaflet.JS

texts:
I have been struggling to use Leaflet (just getting back into coding) and I have done everything I would like to do but add a hyperlink to a pop-up on my map. The current code is:
var marker = L.marker([46.423611, -86.624417], {
        title: "This is a test" ,
      }).addTo(map)
      .bindPopup("<h1> Test Title </h1> <img src=' Link to Image' width='100' height='100' />");

I would like to hyperlink 'Test Title' but each time I try using basic HTML I cannot seem to get it to work. Just a basic hyperlink to a webpage. I would also like to hyperlink the image, but I would stop at just the text for now to save time on this project!


texts: Why doesn't Polygonizer copy line ids to polygonized shapefile?

texts:
I am trying to convert a line shapefile into a polygon shapefile using Polygonizer plugin, and when my new polygons are created all the cells in my id field are "Null", although I have clicked "copy table field from line layer". I have tried both Polygonizer 2 and Polygonizer 2.1.

Any idea why the ids from my line layer are not copied to my polygonized shapefile?


texts: QGIS SAGA Raster Calculator returns incorrect results

texts:
The SAGA Raster Calculator within QGIS produces wrong results probably due to interpolation but works corretly withion SAGA itself. The problem has been discussed here before but sofar I have not found any useful solution.
see also:
   Why does SAGA raster calculator interpolate between cell values?
   Raster calculators in SAGA return incorrect results).
Why does interpolation (?) occur and how can it be avoided? I need the raster calculation to be performed within a QGIS python script.


texts: Leaflet Magnifying Glass and opacity controls

texts:
I have been looking at Leaflet as an alternative to OpenLayers 3 in the event that I need something more lightweight, but I am running into a problem when trying to integrate some of the generously provided plugins. Specifically, I am trying to get Magnifying Glass (https://github.com/bbecquet/Leaflet.MagnifyingGlass/tree/leaflet-0.7.2) and Opacity Slider (https://github.com/lizardtechblog/Leaflet.OpacityControls) to work together. I am trying to use the DeepZoom tilelayer to make an image viewer for large images.

Do these plugins work together? I can get the buttons to show up, but I get t.onAdd errors when I look at the button click event in the console; the map is not altered at all (no creation/removal of the layer for the magnifying glass).

Similarly, I cannot get the slider (jQuery-UI slider) or the buttons from the Opacity Slider to work. It always come back with an error: $(opacity_slider_div) is not a function.

Does anyone know what am I doing wrong?


texts: NAIP Google Earth Engine imageCollection Without Cloud Cover

texts:
I'm trying to gather an imageCollection of an 18 month period of Richmond, California.

Problem 1: extent not working
Problem 2: print() not working

Below is my attempted code that is not working.

// create NAIP image collection within time frame and extent
var dataset = ee.ImageCollection('USDA/NAIP/DOQQ')

.filter(ee.Filter.date('2017-03-01', '2019-06-15'));                .filterBounds(ee.Geometry.Point(-122.35107642463032,37.930370885512346));

// make it so the image is original
var trueColor = dataset.select(['R', 'G', 'B']);
var trueColorVis = {
  min: 0.0,
  max: 255.0,
};
Map.addLayer(trueColor, trueColorVis, 'True Color');

//the information about collection should come up in the console

print(dataset);


texts:
I am working with a shapefile that has many polygons.

How do I add one more field named "area_sqkm" and calculate area for each polygon in the shapefile?


texts: How to increase symbol size of polygons when zoomed out?

texts:
I have many small polygons. How do I increase size of these polygons when zoomed out in QGIS?  As seen in the image below, the symbol appears to be lost when I zoom out.

Notes:
I have created a thematic layer on the polygons.
All layers are vector layers from a database.


texts: Custom color palettes not visible to other layers?

texts:
Full-disclosure: fairly new to QGIS (v3.4), and I am currently styling a map for presentation. I'm trying to match colors within data contained on a layer using "Categorized" styling, so that I can differentiate hierarchy between State/Local/Private/etc roadways, and also do this between discrete layers.

It seems like the way to do this would be to create a "New Palette" of swatches within the "Change Symbol Color" dialogue. Normally, I would expect this would allow me to save a dozen color swatches I like, and then load this palette when styling another item, either within the same layer, or another item on a different layer - but that does not seem to be the way it works. Instead, I find the palette which I've saved, but it is empty with no contents (???). Is this a bug or an environmental variable thing? Do I need to export/save the palette outside of QGIS and load it for each different layer? This seems like a total P.I.A. Anyone else experience this or have a workaraound solution?


texts: How to Make Model with Iterate selected feature run faster?

texts:
This question is related to Fixing Model with Iterate selected feature?

The answer provided by worked great, but it’s really slow. Any way to make it run faster?

My Points Shapefile contain 6707 features and my Flow direction raster is calculated from the SRTM (250m resolution). My data is in GCS in the WGS84 datum.


texts: Applying same color ramp over multiple maps with different values in ArcGIS10?

texts:
I'm doing a rainfall analysis for multiple events in ArcGIS10.1, and I want to apply the same color ramp (with a constant range) for each map. For example, all of the events I have range from zero to roughly 6.5 inches. The color ramp I created is a simple rainbow scheme, with red representing 0 and purpleish representing 6.5. My problem is is that obviously, not all of the events range from 0-6.5. Some go from 0.3-0.7, others from 0.1-1.3, etc. How do I go about applying the 0-6.5 color ramp on each map so that the values are represented with the proper color?

Let me know if this question needs any clarification.

Thanks!

EDIT -- So perhaps I'm doing something wrong.  I made the color ramp as shown (saved it as a layer file and transformed it into a raster) and when I go to import a symbology to one of my interpolated layers  as shown, it stretches the entire color ramp from the lowest value to the highest value, not just the parts of the color ramp it "should". Any advice on what I'm missing?


texts: How to calculate if a turn is required from a set of coordinates?

texts:
I have some polyline data that represents a navigation route generated using Dijkstra's algorithm. Along this route I would like to calculate anytime a turn is required. Ideally, this would be some threshold where the route would require a turn greater than X amount.

Given a polyline, what are some techniques to calculate a turn and it's direction?

I am not using road data from anything like OSM. The data is custom captured coordinates for a traveled path. This means I can't use something that consumes road networks like from OSM data.

I am looking for mathematical techniques to calculate the turns in my custom data. For example, I am currently using two segments of my path. Point1 -> Point2 and then Point2 -> Point3. With these two vectors I am trying to determine if there is an angle in direction change that could be considered a "turn". I am not sure of all the variables I need to take into account here with heading/bearing, etc. Any explanation of what to consider would help.

Current Attempt:

Converting the lat,long data into cartesian coordinate values
Creating two vectors p1->p2 and p2->p3.
Normalizing those vectors
Taking the cross product of those vectors
Cross product > 0 turn right, less than 0 turn left, otherwise straight.


texts: Passing dynamic query parameters to an SQL query string in GeoPandas read_postgis() method

texts:
Up to now, I used to query geospatial tables with psycopg2 and then convert the geometries to my needs using, e.g. Shapely.
Per the psycopg2 documentation, you never format your query string yourself. Instead, you let the .execute() method of your cursor doing the actual job.
import psycopg2
sql_query = ("SELECT * FROM table WHERE id = %(id)s;")

with psycopg2.connect(**params) as conn:
  with conn.cursor() as curs:
    curs.execute(sql_query, {'id': id})
    res = curs.fetchall()

Now, I'd like to use .read_postgis() of GeoPandas to make an equivalent clean and secure call to my database.
But this method is a little bit skinny regarding the query string:

Is there a way to dynamically pass a dictionary of parameters to a call of GeoPandas ' .read_postgis() method? Maybe relying on the same mechanism that psycopg2 is using?
For the moment, my workaround is to take advantage from the psycopg2 way of doing by "mogrifying" the execution:

as such:
import geopandas as gpd
gdf = gpd.read_postgis(
        cursRemote.mogrify(sql_query, {'id': id}), # <- passing this guy instead of my raw query
        connector, geom_col='mygeometrycolumn'
    )

I would be glad to know if it exists a direct GeoPandas' way?


texts: Designing ArcPad Quickforms to Collect Biological Data?

texts:
I am trying to design a Quickform in ArcPad 10.2 to enable efficient capture of biological survey information in the field on a Trimble Juno 3D.

I am currently investigating using the Form Creation Wizard applet guidance used.
I need to be able to pick a large number of plant species within a 'Species' field in the Quickform for each shapefile (each shapefile representing different habitats that I may encounter on a site survey i.e. grassland, woodland, wetland etc.), however, it would appear that only a single Pick List value can be ascribed to each field, in this case my 'Species' field, for each shapefile.

Is there a way that a number of species can be entered into my 'Species' field for each shapefile? There are situations where the habitat may include up to 60 species (i.e. grassland or wetland), although this is quite rare in the type of survey work I currently undertake. Ideally, I would like to have a pick list containing all 7,800 or so plant species in the checklist of British and Irish flora so that I can select from this exhasitive list of all the plant species I could potentially record in each habitat during a site survey. However, I could make do with a much abbreviated list if this would be too onerous and impractical, as the vast majority of these species would be superfluous in most situations.

Also, do you need to be a computer programmer and familar with VBScript to use ArcPad Studio or does it have a WYSIWYG user interface so than non-programmers can use it?

Is there a better way of doing any of this?  I would be very grateful for any help with the above as I have been trying to develop the required 'front end' user interface to run on a Trimble Juno 3D in the field and enable collection of spatially accurate biological data in the field for integration with QGIS back in the office. You would think that there would be an off the shelf solution, but alas it has taken many weeks of research and trial and error to try to develop this with many blind alleys experienced along the way.


texts: ArcGIS Online Groups for Public and Organization Accounts

texts:
I need to share data with public accounts from my organizational account within ArcGIS Online. I don't want to make the data shareable to everyone, just within a group. But when I try to invite the public user to my organization group I get this generic error message:

If I create a group with the public account, I can't invite the organizational user and the organizational user can't find the group.

I may be barking up the wrong tree and this may be a limitation of ArcGIS public accounts, but I've been searching for a while and haven't found anything conclusive to that effect. Perhaps there actually is an error in the back-end?

FYI - I am try to share collector apps in my organization. We have about 7 tablets and I want to have a unique account for each tablet. I was hoping to have public accounts for each tablet. In addition we plan to share certain collector apps with members of the public, once again I don't want to share the app with everyone, but only certain members.


texts:
Is there any way of knowing what Travel Modes are available for directions between two places in Google Maps v3?

Of course, requesting directions for each travel mode and checking which ones return ZERO_RESULTS is a possibility, but it implies four http requests and four google api requests, which I'd very much like to avoid.


texts: Google Earth Engine: I am only interested in the values for a specific day each year, how do I accomplish that?

texts:
I have an imageCollection from CHIRPS, which provides world-wide data over daily precipitation for the past 30 years. I am only interested in the precipitation for the day December 20th each of those 30 years (which should return only 30 points), how would I filter for that day? I have a previously specified boundary as well.

var chirpsColl = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")
                  .filter(ee.Filter.calendarRange(365,'day_of_year'))
                  .filterDate('2010-01-01','2018-12-31')
                  .filterBounds(Dalarna);
//print(chirpsColl);
var chirpsMean = chirpsColl.mean().clip(Dalarna);


texts:
I have a bunch of LAS files that have been delivered with different naming conventions.
I want to rename the files to all have the same uniform naming - based on the easting and northing coordinates of each tiles centroid.
I have generated a text file with the old and new names, and am trying to piece together a python script - but am running into this error

Traceback (most recent call last):
  File "C:\LASTools\Test\ReName_Files.py", line 64, in <module>
    newName = newList[indexOldName]
NameError: name 'indexOldName' is not defined

Here's the script. I am very new to Python and so this could take me hours of research and trial and error (it already has).

import os
folder = "C:\LASTools\Test\\"
lookupTable = open(folder + "LookupTable.txt")
tableList = lookupTable.readlines()
lookupTable.close()

old = "TIFF"
new = "TILE"

newTable = []
for line in tableList:
     newLine = line.strip()
     newerLine = newLine.strip(",")
     newestLine = newerLine.strip("'")
     almostFinalLine = newestLine.replace('"',"")
     finalLine = almostFinalLine.split(",")
     newTable.append(finalLine)

header = newTable[0]
indexOld = header.index(old) - 1
indexNew = header.index(new) - 1

oldList = []
newList = []
for item in newTable:
    if item != newTable[0]:
        valueOld = item[indexOld]
        valueNew = item[indexNew]
        oldList.append(valueOld)
        newList.append(valueNew)

for file in os.listdir(folder):
    nameAndExt = os.path.splitext(file)
    oldName = nameAndExt[0]
    extension = nameAndExt[1]
    print oldName

if oldName in oldList:
    indexOldName = oldList.index(oldName)

newName = newList[indexOldName]

os.rename(folder + oldName + extension, folder + newName + extension)


texts: Is the GPS signal being scrambled near military bases

texts: We are navigating our survey aircraft using standard GPS receiver. Sometimes we fly near the vicinity of US Ramstein military base in Germany and it seems we are losing lock (no solution) in that area when making fast loops or turns. It only happens in this area. Is it possible they are somehow scrambling the signal? I haven't found this documented anywhere...


texts:
I was trying to add a GeoJSON feature implementing leadlet, but with no avail. Cannot figure out what is going on?

<!DOCTYPE html>

    <meta charset="UTF-8">
    <title>Leaflet Demo</title>
    <link rel="stylesheet" href="//cdnjs.cloudflare.com/ajax/libs/leaflet/0.7.3/leaflet.css"/>


texts: Latitude & Longitude values from Simple Voyage Data Recorder [SVDR]

texts:
I am receiving values in the form of 0113.7276,N,10339.0121,E from a Ship Data Recorder. I am not a GIS/ GPS expert and am unable to decipher what do these values mean, because they are out of range for standard latitude and longitude values.

What could these values mean?


texts:
I want to set coordinates of my map in epsg 32749.

var view = new ol.View({
    center:ol.proj.transform([43037,9138591],'EPSG:32749','EPSG:900913'),
    zoom:10.8
    });

var map = new ol.Map({
    target:'map',
    renderer:'canvas',
    view:  view

});

and how to use projs4 to my html?


texts: Convert coordinate system of Lidar or TIN?

texts:
I have Lidar data and I want to create a TIN from it, and then run the TIN Triangle tool to create a polygon shapefile of the triangles.

The only thing is that the Lidar data is in a UTM CS with Meters as the x,y and z coordinates. I want the result shapefile to be in a state plane CS with feet as the x,y and z units.  So I have to convert the CS as some point.

At what point can I do that?  The Project tool doesn't run on either Lidar data or on TINs as far as I can tell. Maybe I could do it as the last step, but does Project work on the Z coordinates of a 3D shapefile?

Thanks!


texts: How to deal with decimal separator when uploading geometries to MSSQL server using shape2sql.exe?

texts:
I'm using Shape2Sql.exe to upload .shp shapefile to MSSQL server 13.0.4502.0 on a Win10 64bit laptop. I've added Shape2Sql.exe.config

<?xml version="1.0" encoding="utf-8"?>
<configuration>
  <runtime>
    <assemblyBinding xmlns="urn:schemas-microsoft-com:asm.v1">
      <dependentAssembly>
        <assemblyIdentity name="Microsoft.SqlServer.Types" publicKeyToken="89845dcd8080cc91" culture="neutral" />
        <bindingRedirect oldVersion="0.0.0.0-14.0.0.0" newVersion="13.0.0.0" />
      </dependentAssembly>
    </assemblyBinding>
  </runtime>
</configuration>

so that I could run the SqlServerSpatial130.dll correctly. Also erased ID Column name field (I found this suggestion on this question).

The geometry uploads correctly but latitudes and longitudes are in this format '524601.70435768971' '6525809.6313226018' and I think is a problem of the decimal separator character on shapefile but I couldnt find where to configure it.

This is the .shp file I'm trying to upload.

Shape2Sql.exe is part of the SqlSpatialTools software packege of www.sharpgis.net.


texts: Inserting an URL link in marker's popup in Folium

texts:
I want to insert a live url link into my marker's popup in Python Folium but I don't know what class and code I need to use.
my code is a simple one:
folium.Marker(location[i], popup="www.dingsbums.lu").add_to(marker_cluster)

I just want the url (here 'www.dingsbums.lu') to be a link, not string...


texts: Create a function that will convert CSV to other types of files

texts:
I'm writing a function with my schoolmates that will convert various files (like shapefile to CSV, CSV to TopoJSON etc, we focus on geospatial data) - we are trying to learn Python. My part right now is to write two functions for CSV format (to read and write file) - we write two for each format and then we use them in the main function. This is what I came up with for CSV:

import pandas as pd
import geopandas as gpd

#read file
def from_csv(file, x, y):
    df = pd.read_csv(file)
    gdf = gpd.GeoDataFrame(df, geometry = gpd.points_from_xy(df[x], df[y]))
    return gdf

#write file
def to_csv(gdf, filename):
    gdf.to_csv(filename + '.csv')

I tried to use Pandas to read file and then used GeoPandas to create gdf from df - I want user to write which columns are longitude and latitude (to set it as a geometry while creating the GeoDataFrame). Also I know that DataFrame may have coordinates in WKT format. I don't really know how to even start with including that.

To give you an idea - here is a little snippet of our function, I cut out a lot and only left CSV, trying to make it as short as I can. There are more formats that we already included - like GeoJSON. All other functions work (like converting GeoJSON to shapefile, no problems with that). I tried to convert GeoJSON to CSV and it worked, but I cannot convert CSV files to any other format.

valid_formats = ["csv"]
valid_extensions = [".csv"]

def convert(file, target_format, x = None, y = None):
    filename, file_extension = os.path.splitext(file)
    if file_extension in valid_extensions and target_format in valid_formats:
        if file_extension == ".csv":
            if x != None and y != None:
                gdf = from_csv(file, x, y)
        if target_format == "csv":
            to_csv(gdf, filename)
    else:
        print("wrong input data")


texts: I installed QGIS 3 as instructed and the Vector processing tab is missing. I did uninstall all previous versions and Python.


texts: Securing WMS against unauthorized access?

texts:
In a project the customer wants the images of the WMS to be secured against unauthorized access.

What are the common solutions to realize that?

Are they compatible with the OGC-standards?


texts:
Is there any possibility to make a SD file without using ArcGIS Desktop?

I have to upload my service on a service manager and I am using PostgreSQL as a database.


texts: QGIS: How do I write a CASE expression (field calculator) that doen't calculate/change the existing field information if not true?

texts:
I'm using a CASE WHEN expression in the field calculator on an existing field which already contains information. The expression is comparing two fields. If these fields match, the result should be 0. It looks something like this CASE WHEN "TD_sum" = "Qup THEN" '0' END. The field I'm editing already contains values which should stay the same if the CASE isn't true. So how do I alter the expression so it only calculates the matching fields? Is it possible to add ELSE? What do i write as the output?
I know that I can simply select the matching fields with select features and then calculate them, but maybe this calculation will be necessary for a model and I want as little individual steps as possible.


texts: Creating relationship class for many to many relationship

texts:
In Arcmap 10.7, I am trying to set up a relationship class that will link parcels to owners. The end goal is to be able select a parcel and see the owner, or multiple owners associated with it.

Some parcels have multiple owners and some owners have multiple parcels. I assume this means I need a many to many relationship. I have a feature class of parcel polygons, and a stand alone table with owner information. They are both in the same geodatabase. I have tried using parcels as origin, and also using owners table as origin. I am using a simple relationship both ways. I choose no message propagation since I am unsure what does (how edits are passed along?).

I can never get any records to show up in the related tables. I believe the keys I am using to link them are valid, both tables have many overlapping values. Is there a reason I can't get any related records to show up? This was asked before in Records in related table do not appear when using the 'Related Tables' button in attribute table post but was never answered.

I saw other posts describing using overlapping polygons to solve the multiple owners problem, but that isn't what I would like to do if possible. Would like to just have all the owners associated to the one parcel shape in the feature class.


texts:
Script is simple made for this case. Creating new shapes based on known arguments.

My code:

import arcpy

#variables
fc_1 = 'fc_1'
fc_2 = 'fc_2'
fc_3 = 'fc_3'
fc_4 =  'fc_4'
war_1 = "highway = 'path'"
war_2 = "highway = 'motorway'"
war_3 = "natural = 'water'",
war_4 = "leisure = 'pitch'",
suffix_1 = '1'
suffix_2 = '2'
suffix_3 = '3'
suffix_4 = '4'
temp = 'temp_layer'

def select_save(fc, temp, war, suffix):
    arcpy.MakeFeatureLayer_management(fc, temp)
    arcpy.SelectLayerByAttribute_management(temp, "NEW_SELECTION", war)
    arcpy.CopyFeatures_management(temp, '{0}_{1}'.format(fc, suffix)

arcpy.env.workspace = r'gdb_path'

select_save(fc_1, temp, war_1, suffix_1)
select_save(fc_2, temp, war_2, suffix_2)
select_save(fc_3, temp, war_3, suffix_3)
select_save(fc_4, temp, war_4, suffix_4)

So how to get rid of those copy-pasted function calls?
In my current project I need to call my function 20-30 times. It looks horrible to see this wall of text.


texts: How to change GeoServer memory limits for Cross-layer filter in Linux?


texts: Reprojecting from Lambert Conformal Conic into EPSG:3857?

texts:
I'm trying to reproject a plain png image file into EPSG:3857 for use in an OpenLayers application. I'm pulling the image from the NWS site at: http://rapidrefresh.noaa.gov/HRRR/for_web/hrrr_ncep_jet/2015032517/full/cref_sfc_f00.png

According to that site the map parameters are:

    :CEN_LAT = 38.5f ;
    :CEN_LON = -97.5f ;
    :TRUELAT1 = 38.5f ;
    :TRUELAT2 = 38.5f ;
    :MOAD_CEN_LAT = 38.5f ;
    :STAND_LON = -97.5f ;
    :POLE_LAT = 90.f ;
    :POLE_LON = 0.f ;

Below are corner points for WRF-HRRR mass points, SW-NW-NE-SE

    :corner_lats = 21.13812f, 47.84364f, 47.84364f, 21.13812f;
    :corner_lons = -122.7195f, -134.0986f, -60.90137f, -72.28046f;

I've referenced Making geotiff and reprojecting to EPSG:4326 using GDAL? . I've cropped the frame down to just the mapped area and I've tried running cs2cs to get the map units I need to  run gdal_translate and gdalwarp but I can't seem to get the correct (WebMercator) image out of gdalwarp. I think my problem is that I don't have the projections quite right in cs2cs. I'm running:

 cs2cs +init=epsg:3857 +to +proj=lcc +lat_1=38.5 +lat_2=38.5 +lat_0=38.5 +lon_0=-97.5 +k_0=1 +a=6371200 +b=6371200 +units=m +no_defs
 -134.0986 47.84364
 10996563.15    1837306.24 0.00
 -72.28046 21.13812
 10996629.05    1837356.50 0.00

Then

 gdal_translate -a_srs "+proj=lcc +lat_1=38.5 +lat_2=38.5 +lat_0=38.5 +lon_0=-97.5 +x_0=0 +y_0=0 +k=1.0 +a=6378137 +b=6378137 +units=m +no_defs" -a_ullr 10996563.15    1837306.24 10996629.05  1837356.50 cropped2.tif cropped2_translated6.tif

Then finally:

 gdalwarp -t_srs "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=-90.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs" cropped2_translated6.tif cropped2_warped.tif

My final warped image comes out with N at lower left S at Upper Right E at Upper left and W and lower right so the US looks like it is falling over backwards.

A gdalinfo of this warped file shows:

 gdalinfo cropped2_warped.tif
 Driver: GTiff/GeoTIFF
 Files: cropped2_warped.tif
 Size is 850, 917
 Coordinate System is:
 PROJCS["unnamed",
  GEOGCS["unnamed ellipse",
    DATUM["unknown",
        SPHEROID["unnamed",6378137,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]],
 PROJECTION["Mercator_1SP"],
 PARAMETER["central_meridian",-90],
 PARAMETER["scale_factor",1],
 PARAMETER["false_easting",0],
 PARAMETER["false_northing",0],
 UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]
 Origin = (10012443.303476206958294,7747.939850974368710)
 Pixel Size = (0.072684385553582,-0.072684385553582)
 Metadata:
   AREA_OR_POINT=Area
 Image Structure Metadata:
   INTERLEAVE=PIXEL
 Corner Coordinates:
 Upper Left  (10012443.303,    7747.940) (  0d 3'24.09"W,  0d 4'10.56"N)
 Lower Left  (10012443.303,    7681.288) (  0d 3'24.09"W,  0d 4' 8.41"N)
 Upper Right (10012505.085,    7747.940) (  0d 3'22.09"W,  0d 4'10.56"N)
 Lower Right (10012505.085,    7681.288) (  0d 3'22.09"W,  0d 4' 8.41"N)
 Center      (10012474.194,    7714.614) (  0d 3'23.09"W,  0d 4' 9.49"N)
 Band 1 Block=850x3 Type=Byte, ColorInterp=Red
 Band 2 Block=850x3 Type=Byte, ColorInterp=Green
 Band 3 Block=850x3 Type=Byte, ColorInterp=Blue

What am I doing wrong?


texts: Using URL of JSON file on local drive with OpenLayers support?

texts:
When i save json file in web apps then using http://localhost:8080/WGSClient/js/object.json this path file gets open.

But when i keep it outside web apps like local folder then using same url format json file is not accessible.

It showing error that OpenLayers doesn't support this url format.

So what is the proper url which OpenLayers support or i should keep my json file in web apps?

What would be the URL of a JSON file on the local drive?


texts: HRESULT Error 0x80041018 when using IEnumDataset.Next()

texts:
Similar to HRESULT E_FAIL Error -2147467259 when using IEnumDataset.Next(), I'm writing a Desktop Add-In that scans our directories and pulls in georeferenced raster images that contain a point geometry. Whenever I get to the image shown below, I get the HRESULT error 0x80041018 (or -2147217384).

I can pull this image into ArcMap and it renders just fine so I know it isn't corrupt. I tried removing the image from the folder and seeing if it would run successfully but it threw the same error on one about 20 images previous to it (even though it had previously had no problem).

When I looked up the error code here it says:

  FDO_E_INVALID_CONNECTION_POINT_GEOMETRY     -2147217384     The geometry for a complex junction connection point is invalid.

Not quite sure what to make of this error.

Code:

//how RWSF and IPoint are declared in other methods
Type factoryType = Type.GetTypeFromProgID("esriDataSourcesRaster.RasterWorkspaceFactory");
IWorkspaceFactory RWSF = Activator.CreateInstance(factoryType) as IWorkspaceFactory;

IPoint point = new PointClass();
point.SpatialReference = NAD83;
point.PutCoords(longitude, latitude);

//method that uses them
private static void PullImages(IPoint centerpoint, string directoryPath)
{
    IWorkspace ws; IDataset dataset; IRasterLayer rasterLayer;
    ws = RWSF.OpenFromFile(directoryPath, 0);
    IEnumDataset enumDataset = ws.get_Datasets(esriDatasetType.esriDTRasterDataset);
    dataset = enumDataset.Next();

    while (dataset != null)
    {
        IEnvelope rasterEnvelope = (dataset as IGeoDataset).Extent;
        rasterEnvelope.Project(NAD83);
        if ((rasterEnvelope as IRelationalOperator).Contains(centerpoint))
        {
            rasterLayer = new RasterLayerClass();
            rasterLayer.CreateFromDataset(dataset as IRasterDataset);
            MXD.AddLayer(rasterLayer);
        }
        dataset = enumDataset.Next();
    }
}


texts:
My aim is to plot a rasterstack as a single image using ggplot2. First, I apply ggplot on a single layer raster using the following code:
#Load packages
library(raster)
library(ggplot2)

#Create a random raster layer
set.seed(11)
r <- raster(nrows = 10, ncols = 10, res = 30, xmn = 267195, xmx = 267375, ymn = 4016985, ymx = 4017135)
r <- setValues(r, runif(ncell(r), min = -10000, max = 10000))
crs(r) <- "+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
plot(r) #Produce simple plot

rdf <- as.data.frame(r, xy=TRUE) #Convert raster to data.frame
names(rdf)[3] <- 'magnitude' #Name value column
head(rdf)

Now I apply gglpot()
ggplot(data = rdf)+
  geom_raster(mapping=aes(x=x, y=y, fill=magnitude))+
  scale_fill_gradientn(colours= rev(terrain.colors(10)), name='Magnitude')

It successfully prints. Now I create a raster stack with 10 layers and call this ndvi.
rr <- lapply(1:10, function(i) setValues(r, runif(ncell(r), min = -10000, max = 10000)))
ndvi <- stack(rr)
crs(ndvi) <- "+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
plot(ndvi)

How could I plot this stack using ggplot2? Do I create a data.frame but then how do I supply fill argument in aes() and facet_wrap()?


texts: Mapquest directions API in ArcGIS javascript webmap application

texts:
I'm trying to build a webmap application using the ArcGIS JavaScript API.  A functionality I would like to include is driving directions from a current GPS location to a feature in my feature service.

I don't want to use Esri's ArcGIS Online service because it consumes credits.  I created a geoprocessing routing service which is hosted in ArcGIS for Server, but getting it to work has a whole other set of complications.  My last resort is to try and use a driving directions API like that provided by Mapquest.

Is there a way to incorporate this type of API in an ArcGIS JavaScript API, and if so, how?


texts: Can Quantum GIS replace Proprietary GIS?

texts:
I’m an MSc GIS student in the final year of my studies. I’m researching for my dissertation project which is entitled ‘Can Quantum GIS replace Proprietary GIS in large organisations?’

I’m primarily seeking the views of GIS professionals/Managers who have taken the decision to replace Proprietary GIS (PGIS) with Open Source GIS OSGIS (Not necessarily Quantum GIS) or have decided to stick with PGIS. I’m also interested in the views of GIS users who have considered migrating to OSGIS. The main questions I have are as follows:

•   Have you considered switching to OSGIS?

•   Why did you decide to migrate to OSGIS or to retain PGIS?

•   If you migrated to OSGIS what issues did you encounter?
-   Are you still encountering issues, if so what?
-   With your experience in mind, what are the pro’s and con’s of OSGIS
-   Has your GIS functionality/capability increased as a result of migrating to OSGIS?

•   If you decided to retain PGIS, what was your reasoning behind this?
-   What were your concerns relating to OSGIS
-   With your experience in mind, what are the pro’s and con’s of PGIS
-   What improvements, if any, would make you consider migrating to OSGIS?

Feedback from all GIS users is most welcome. I have chosen not to produce a survey as I want the responses to be open ended, allowing users to add as little or as much feedback as they wish. I’m also hoping users will raise issues that I previously hadn’t encountered.

Thanks in advance.


texts: Counting within buffers using ArcGIS Desktop?

texts:
I am working on a project where I am trying to map affordable housing and bus stops.

I have used the buffer tool to find out how many bus stops are there within a one mile radius of each affordable house but is there a way to count the bus stops in that one mile radius?

I want to say near X affordable house, there are Y bus stops. Right now what I have is a circle indicating the one mile radius but no count of the number of bus stops in that area.


texts: How do I add a lot of data to a shapefile in QGIS?

texts:
I have some nice shapefiles from Agesic in Uruguay, I need to add to them some information (not as shapefile) that I got from Dinamige to the attribute table. These institution are not the point, the point is that I have a lot of new information to add to the shapefiles.

I was able to add new columns, they look very nice and empty, so I want to copy and paste the new information there and not write it one by one. I have just one attribute that is the same, all the other stuff is new information. The thing is I'm a new QGIS user and I don't understand how Python Console works, that is like out of my league. I was searching someone with the same problem but I don´t find it. Does anyone know how to do this?


texts:
I am trying to use a custom button in Openlayers but is does not work. what am i doing wrong?

var custom_button = new Openlayers.Control.Button({
            displayClass: 'olControlCustomButton',
            trigger: codigo2,
        })

        var control_panel = new Openlayers.Control.Panel({});

        control_panel.addControls([
            custom_button
        ])

        map.addControl(control_panel);


texts: How can I generate a range of pantones for a map?

texts:
When I am given a map (or make one...) I often like to get an idea of the range of colours present. One thing this does is it allows me to quickly visualize contrast (helps build maps for colour blind readers).

I was wondering what some options are to automatically or easily generate a range of pantones for any given map? There may be a different process when dealing with vector based maps, like a PDF , PS/EPS, or a SVG.

I wonder if there is a tool or plugin for GIMP that will scan an image, get the rgb/cmyk value for every pixel and summarize with a count showing the occurrence...

An opensource suggestion would be ideal. I use GIMP a fair bit, but only for literal design/layout and have not delved into image analysis.

I have ArcGIS 10 but without spatial analyst so I have no raster tools available in that package.

Sorry if this is on the wrong SE site, maybe better for photo.SE?


texts: Combine overlapping rasters and include cells with no data

texts:
I'd like to combine about a dozen rasters into one raster that contains the "summed" cell values where overlaps occur.  In doing so, I need to be able to retain cells that either have no data, or that do not overlap with any other rasters.

FYI, Cell values are stored in an integer field in each raster.

I have tried the Raster Calculator, however it doesn't seem to do well with the "No Data" cells- it only creates summed cells where rasters overlap.

I have also tried the "Cell Statistics" tool, however the output feature it generates doesn't have any data.


texts: ESRI License Initialize & Binding Failure

texts:
Despite following guidance about how to first bind, then initialize a licence for a standalone WinForm uses some arcobjects, I still get the "80040111 ClassFactory cannot supply requested class" error once I call GxDialogClass().  When debugging, the status object beneath comes back as licence checked out - so it successfully binds AND gains a licence.  Yet the 80040111 indicates that I haven't.

We are running ArcGIS 10.2.1.

ESRI.ArcGIS.RuntimeManager.Bind(ProductCode.Desktop);
IAoInitialize aoInitialize = new AoInitializeClass();
var status = aoInitialize.Initialize(esriLicenseProductCode.esriLicenseProductCodeAdvanced);
IGxObjectFilterCollection dialog = new GxDialogClass();

Any ideas?

Matthew


texts: ModelBuilder make no selection a valid selection?

texts:
Using ArcMap 10.3.1

What I am trying to achieve using ModelBuilder...

I want end users to select either 'roads', 'protected lines' or a combination of the two and then have the selected features merged into a new feature.

The problem I am having is if I select 'roads' but no 'protected lines' (or vice versa) the model assumes no selection means use the entire data set when I in fact want it to mean no selection/no data to merge.

I've tried incorporating the 'get count' tool but the result is the same as above (maybe I'm using it wrong?)

Image below of my model.


texts: Adding sigmoidal curve function in ArcGIS raster calculator?

texts:
I would like to add a sigmoidal function to create a factor map In ArcMap. The function has to decrease with an increase in values.

For example: The suitability should be 1 for a distance of 250m and then it should decrease as a cosine function until 1250m.

250m has highest suitability value 1 and it decreases to 0 for distances more than 1250m

Could someone explain how to do this using raster calculator?


texts: Editing already published web mapping applications

texts:
Is there a way to update the web mapping application on ArcGIS Online?

I have an application that I want to add layers to from our REST services and I can't seem to find a way to add it. It does not give me the option to edit the application. It gives me the options to view item details and view application.

I do not have to original web map for the application - so I wouldn't be able to republish it.  It was originally published in 2015.


texts:
This is perhaps a silly question but can ArcGIS Explorer open a .mxd file?

I can't seem to find a definite answer online and I can't download it until I know the answer.


texts: QGIS: Cut line feature by the nearest point

texts:
I got a line shapefile and a point shapefile (or a csv file containing coordinates of the point), and I want to clip this line feature by those point in the 'nearest' way.

So far, I have always used edit function in QGIS GUI to clip line feature manually:

I am wondering are there QGIS built-in functions could handle this problem (for example, input line feature and point feature, output cut line feature (new file is even better).)? or could I cut line feature in this way using PYQGIS command or user-written python function (even better)?

Update#1

I used PyQGIS code from @xunilk to split line feature showed in figure above, and named my input line and point shapefile as 'my_border' and 'my_point' respectively, but the results looked like:

The result is not what I expect.

Besides, I followed the way in the answer provided by @Micha in splitting a lines layer using points, trying to split line by GRASS.

Because my point is not on the line, I have to run the following module to get 'error point':

v.distance
v.patch
v.clean

I think there are no problem use the first 2 module, but after I used v.clean:

v.clean --overwrite input=tcall@nl output=tcerrpoint tool=break

I just got a line feature, not error points:

Were there something I missed or did wrong in the procedure?

Update#2

I used the updated code provided by @xunilk, which trying to handle multiple cut points. I still used my line and point shapefile as input layers, and the following figures are their attribute tables.

The output is:

Still not what I expected, I am wondering that maybe it outputted strange line feature because the input line has too many segments.

If the PyQGIS script being modified to capable of handling selected feature of input line shapefile, and before run the script we select the nearest line segment to one point, maybe the output can be right?

Update#3

This time I edited input line shapefile 'my_border' to a line with only one attribute:

and used PyQGIS code in the edited answer of @xunilk:

registry = QgsMapLayerRegistry.instance()

my_border = registry.mapLayersByName('my_border')
my_point = registry.mapLayersByName('my_point')

feats_my_border = [ feat for feat in my_border[0].getFeatures() ]
feats_my_points = [ feat for feat in my_point[0].getFeatures() ]

points = [ feat.geometry().asPoint() for feat in feats_my_points ]

my_tuplas = []

for point in points:
    tupla = feats_my_border[0].geometry().closestSegmentWithContext(point)
    my_tuplas.append(tupla)

my_line = [ feat.geometry().asPolyline() for feat in feats_my_border ]

new_line = []

k = 0

try:
    for i, line in enumerate(my_line):
        for j, point in enumerate(line):
            new_line.append(point)
            if j == my_tuplas[k][2] - 1:
                new_line.append(my_tuplas[k][1])
                k += 1

except IndexError:
    pass

k = 0

lines = [ [] for i in range(len(points) + 1) ]

epsg = my_border[0].crs().postgisSrid()

uri = "LineString?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'split_line',
                           'memory')

prov = mem_layer.dataProvider()

try:
    for point in new_line:
        lines[k].append(point)
        if point == my_tuplas[k][1]:
            k += 1
            lines[k].append(point)

except IndexError:
    pass

for i, line in enumerate(lines):

    feat = QgsFeature()
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPolyline(line))
    prov.addFeatures([feat])

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

This time the attribute of output 'split_line' looked rational because I inputted 8 point to cut line:

But no matter what I tried including set layer CRS, zoom to layer, show in overview, etc., I could not see this output file. Besides, I saved this output file as a new ESRI shapefile and checked its size, which was 208 kb (input 'my_border.shp' is 21 kb), all seemed to be rational, but I just can't see the output.


texts: Cannot log in with different ArcGIS Portal user account using ArcGIS Python API

texts:
I have an ArcGIS Portal log in account (viewer account) in an enterprise GIS infrastructure and want to use the ArcGIS Python API to connect with a separate admin account to be able to produce reports on the ArcGIS Portal instance (layers being used in web maps, user last login stats etc).

When trying to connect using an admin account, it automatically connects with my viewer account (I am assuming my account is somehow automatically authorised within my network).

How am I able to explicitly log in with a different account using the ArcGIS Python API?

My script is below:

from arcgis.gis import GIS
import getpass

password = getpass.getpass("Enter password: ")
gis = GIS("https://myserver", "myadminaccount", password)
print("Connected as: " + gis.properties.user.username + " on " + gis.properties.customBaseUrl)
print(gis.properties.user)


texts: How to use fanout to accomplish batch clip in FME?

texts:
Following up on this question. I am trying to implement the fanout method to implement a batch clip that was mentioned as an option.

I am having issues with trying to batch clip in FME. This is my workspace:

I feed the Clipper 3 shapefiles (each has only one part) to clip the larger clippee shapefile  (that has 59 parts). How do I set up the fanout in the writer so that I write 3 separate shp files? Currently, it writes one shapefile with 15 parts.

What I tried was to expose the fme_basename of "All" ( the 3 shapefiles read in as a single merged feature type), and then fanout by fme_basename name in the writer.


texts: How to return count of popups opened up by clicking on a mapbox gl map feature?

texts: I’ve a point dataset with overlapping points. For each feature a popup is created automatically. When I click on a point with other overlapping point on it all popups are opened up and it’s messed up. Is there a way to count overlapped points to create some kind of a list for those cases?


texts: How can I change region style in MapBasic

texts:
I've problem with my code in MapBasic for change region style in MapBasic.

I've created voronoi and I want to select that voronoi and change the style of it to yellow border and no fill. I already try a code to change that color which I got it from other source but it takes a long time more than I change it manually from MapInfo and when I try to run it becomes not responding.

This is my code

Sub Voronoi

    Create Table "voronoi"
        (Block_no Char(15),Remark Char(10),Type_Palm Char(10),Ha Float)
        File "D:\voronoi.TAB"
        TYPE NATIVE Charset "WindowsLatin1"

    Create Map For voronoi CoordSys Earth Projection 1, 104

    Create Index On voronoi (Block_no)

    Add Map Layer voronoi

    Set Map Layer 1 Editable On

    Create Object As Voronoi from sensus Into Table voronoi

    Update voronoi Set Ha = Area(obj, "sq m")

    select * from voronoi
        where Ha >= Val(inputnumber)
        into Selection

    Call ChangeColor

    Browse * from Selection

End Sub

Sub ChangeColor

    Fetch first from Selection
    x = 1
    Do
        oUpdateObj = Selection.obj

        b1 = Makebrush (1, 0, 16777215) '<--- Change values for your style
        p1 = MakePen(1, 2, 16776960) '<--- Change values for your style
        Alter Object oUpdateObj Info OBJ_INFO_BRUSH, b1
        Alter Object oUpdateObj Info OBJ_INFO_PEN, p1

        Update Selection
            Set Obj = oUpdateObj
            where RowID = x

        x = x + 1

        Fetch next from Selection
    Loop while x <= TableInfo(Selection,TAB_INFO_NROWS)

End Sub


texts: Transform coordinate from rotated lat/lon to normal lat/lon WGS84

texts:
I have access to an OPeNDAP service that allows me to peek into a humongous ocean wave dataset. It takes latitude and longitude as inputs (in fact even a range of the two) - but the coordinate system is a rotated one (presumably because the dataset focus is the North Sea and the Arctic):

grid_mapping_name: rotated_latitude_longitude
grid_north_pole_longitude: 140.0
grid_north_pole_latitude: 22.0
earth_radius: 6371000.0
proj4: +proj=ob_tran +o_proj=longlat +lon_0=-40 +o_lat_p=22 +R=6.371e+06 +no_defs

I want to query specific coordinates, but I do not know how to translate these coordinates (i.e. WGS84 lat/lon from Google Maps) into coordinates for this service. Anyone have an idea how to do this?

For example, let's say I want to transform the coordinate (57.93, 3.63) (somewhere out in the North Sea).

Any solution is welcome, using pyproj or something completely different.

EDIT:
The OPeNDAP dataset can be found Here - it sadly does not offer a whole lot of information, but snippets of data can be extracted as ASCII.

EDIT2:
This is the code I've tried, but it seems the proj4 string is not correctly interpreted:

from osgeo import osr
src = osr.SpatialReference()
tgt = osr.SpatialReference()
src.ImportFromEPSG(4326)
tgt.ImportFromProj4('+proj=ob_tran +o_proj=longlat +lon_0=-40 +o_lat_p=22 +R=6.371e+06 +no_defs')

transform = osr.CoordinateTransformation(src, tgt)
coords = transform.TransformPoint(-122, 46)
x,y = coords[0:2]

print(coords)

Returns:

(-13580977.876779377, 5780349.220256355, 0.0)
ERROR 6: No translation for an empty SRS to PROJ.4 format is known.
Traceback (most recent call last):
  File "coordtrans.py", line 24, in <module>
    coords = transform.TransformPoint(-122, 46)
  File "C:\tools\Anaconda3\envs\testEnv\lib\site-packages\osgeo\osr.py", line 1283, in TransformPoint
    return _osr.CoordinateTransformation_TransformPoint(self, *args)
NotImplementedError: Wrong number or type of arguments for overloaded function 'CoordinateTransformation_TransformPoint'.
  Possible C/C++ prototypes are:
    OSRCoordinateTransformationShadow::TransformPoint(double [3])
    OSRCoordinateTransformationShadow::TransformPoint(double [3],double,double,double)

Run in a conda environment with gdal from conda.


texts: Getting index number of list item using Python?

texts:
this gives me a list of all 55 layers in my mxd

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "Layer")[0]
for lyr in arcpy.mapping.ListLayers(mxd, "*", df):
    print lyr.name

What I want to do is get the index number of a specific layer.  For example, I want to look for the layer named 'Water' and I want the script to tell me "Water' is at position 50, so I know to look for it at the bottom on the TOC.

How can I do this?


texts: Producing building shadows using ArcGIS Desktop?

texts:
How would you generate building shadows from building footprints?
I have attribute values for the heights of buildings.
So would preferably like to have the shadows' sizes vary based on that.
The best suggestion I've found using ArcGIS is drop shadows: http://blogs.esri.com/Support/blogs/mappingcenter/archive/2008/10/21/3d-effect-for-a-building-footprints.aspx


texts: How can I apply multiple expressions to a field in QGIS?

texts:
I'm working expanding abbreviated names in a shapefile. I have written the following expansions, and would like to learn how to write a field expression that would do many expansions in one.

replace( "DNRNAME",'SF','State Forest')
replace( "DNRNAME",'FMA','Fish Management Area')
replace( "DNRNAME",'FT','Fire Tower')
replace( "DNRNAME",'SP','State Park')
replace( "DNRNAME",'CFL','Chesapeake Forest Lands')
replace( "DNRNAME",'NRMA','Natural Resource Management Area')
replace( "DNRNAME",'NEA','Natural Environment Area')
replace( "DNRNAME",'WMA','Wildlife Management Area')
replace( "DNRNAME",'HCF','Heritage Conservation Fund Site')

How can I apply all of these replace conditionals at one time?


texts: Why do my geowebcache tiles not work in openlayers client?

texts:
var options = {
        resrictedExtent: restrictedBounds,
        maxExtent: restrictedBounds,
        units: "m",
        projection: "EPSG:31467",
       resolutions:[305.74872463815512,152.87436231907702,76.437181159538511,38.218590579769391,19.109295289884642,9.554647644942321,4.7773238224711605, 2.3886619112355802, 1.1943309556178034, 0.59716547780890172],
       maxResolution: "auto",
        controls: []
    };

map = new OpenLayers.Map( 'map', options );

These are the map options which work with geoserver wms fine. I want to serve geowebcache layers using geoserver 2.5 with the integrated GWC. I have created a cached layer using the resolutions as they are defined above. In GWC the layer preview works OK and displays the layer as expected, but when I integrate the layer into the openlayers client, pink tiles are displayed. The resolutions seem indeed to differ from those defined by the normal wms layers. Have I forgotten something? Why don´t the resolutions match? The DPI is 96 by the way.

The URL of a wms tile which displays fine in zoomlevel 8 is:

SRS=EPSG%3A31467&BBOX=3542676.9162395,5865437.4059784,3563009.2064279,5885769.6961669

If I change the URL to point to the GWC layer but keep the BBOX extents the same I get the following error:

400: Unable to match requested SRS EPSG:31467 to those supported by layer


texts: Clustering without aggregation with QGIS Point Displacement

texts:
I've seen the following map on twitter and am looking to create the same effect for point displacement.

I'm familiar with grid and circle point displacement as well as with clusters point to aggregate points into one (depending on the scale) but I'm looking for something in the middle.

  clustering point without aggregation, like bunch of grapes

I am tending to achieve this result:

How can I perform that in QGIS?


texts: Some gridlines missing in QGIS Print Composer

texts:
I'm very new to this so please excuse my ignorance.

I've been working with OS data and my own shape files and had no problem with print composer.  Suddenly my OS gridlines are not appearing properly.  Some are there, some aren't.  The map is fine until I go into print composer.

What have I done wrong?

Thanks


texts: How to add a KML vector from a string variable

texts:
I have a variable that contains KML data how so I tried to do some research and come up with this code

var kmlString = result;

var features = new ol.format.KML().readFeatures(kmlString, {
    dataProjection: 'EPSG:4326'
    });

var KMLvectorSource = new ol.source.Vector({});

var KMLvector = new ol.layer.Vector({
    source: KMLvectorSource,
    visible: true
    });

KMLvectorSource.addFeature(features);
map.addLayer(KMLvector);

but I don't know why it doesn't work, I get error message: a.addEventListener is not a function

My KML string is:

<Polygon><extrude>0</extrude><tessellate>0</tessellate><altitudeMode>relativeToGround</altitudeMode><outerBoundaryIs><LinearRing><coordinates> list of coordinates </coordinates></LinearRing></outerBoundaryIs></Polygon>


texts: What to do in ModelBuilder after adding layers and dragging tool?

texts:
I'm trying to create a model that would Merge two layers and then do a Summarize to a common attribute/field between them.

The process I'm doing is getting into the ModelBuilder then adding the 2 layers, then Drag the merge tool into it, then I'm lost!

Is there any possible way to correct my steps or add any steps?


texts: How to connect to (& edit) PostgreSQL database using ArcGIS Desktop without ArcSDE?

texts:
After trying in vein to get spatialkit and zgis working with my postgresql database in arcgis 9.3 and 10 , I was wondering if anyone has had any luck with these products.

My postgresql database has been serving geospatial data for 2 years through geoserver and can be edited in quantum gis fine....so why not in ArcGIS?

Does anyone have any other recipes for using postgresql data with arcgis but WITHOUT ArcSDE?


texts: Sync ArcMap, AGOL, Collector without server

texts: Is there a way to sync one file gdb across catalog or arcmap, agol, and collector without using a server? What is the best way to track changes to a file gdb when people are editing it in the filed (collector) and remotely (desktop) and it is stored in catalog as well as online?


texts: Adding PDF to backside printing for Data Driven Pages?

texts:
Is there an easy way to add a PDF to the backside of each map printed using data driven pages?

I know I could print one sided and then reinsert the printed pages upside down, and print the PDF on the reverse side. However, we may deliver this to our client digitally and would like it to be setup for them.

For example, I have data driven pages set up for 50 parcels. Each map printed will show a single parcel's data. I need a generic PDF to print on the backside of each of these 50 maps.


texts:
Is there a way to format the way attribute values are displayed in the attribute table?

For example, I would rather display the following columns to show comma-separated thousands and without the zeroes after the decimal point.

So, instead of "2300000.00000", it would display as "2,300,000" instead.


texts: Using WFS layers in Virtual Layer crashing/too slow QGIS

texts:
I am trying to create a virtual layer with the following code, where the layers are provided by a WFS:

select geometry from (
select geometry from "Ancient Woodlands (England)" union all
select geometry from "Areas of Outstanding Natural Beauty (England)" union all
select geometry from "Country Parks" union all
select geometry from "Flood Map for Planning (Rivers and Sea) - Flood Zone 3" union all
select geometry from "Flood Map for Planning (Rivers and Sea) - Flood Zone 2" union all
select geometry from "Local Nature Reserves (England)" union all
select geometry from "National Parks (England)" union all
select geometry from "RAMSAR (England)" union all
select geometry from "Sites of Special Scientific Interest (England)" union all
select geometry from "Special Areas of Conservation (England)" union all
select geometry from "Special Protection Areas (England)" union all
select geometry from "SSSI Unit (England)" union all
select geometry from "World Heritage Sites" union all
select geometry from "Registered Parks and Gardens"  union all
select geometry from "Registered Battlefields")
as Constraint_Area

The WFS layers are set to load only the visible extents (restrictToRequestBBOX='1'). The layers load quickly individually but the virtual layer takes forever causing QGIS to stop responding. Is there a way to perform this query successfully?


texts: How to debug settings page of custom widgets?

texts:
I'm developing a setting page for my widget. To debug it, I consulted the guide again and found this advice:

  During development, use http://your host/webappbuilder/?appid=stemapp
  to directly access your widget. This URL helps test the configuration
  page of the widget.

So I tried to open my builder with ?appid=stemapp query string. I got the list of my apps as though there was no query string. The only difference is that when I open an app from this page the loading lasts forever, not few seconds as when opened from the normal WAB root. This doesn't seem to be the direct access for my widget. Google on web appbuilder "appid=stemapp" found only the guide and two pages citing it, so it's possible that it is just a documentation error and the direct widget access either doesn't exist or is somewhere else.

So, do I have any better options for setting page debugging than app edit > widgets? While editing the app, especially early stages of the debugging are frustrating, Create widget setting page error:widget/MyWidget/Widget is not very helpful.

I reposted this on the ESRI support forum.


texts: Export planet OSM file from Rails port using Osmosis

texts:
In order to create a local OSM server and connecting JOSM to this, we've installed and configured Rails port openstreetmap/openstreetmap-website. The next step was to insert a sample OSM file into Rails' Postgresql databae, using osmosis via the following command:
osmosis --read-xml file="test.osm" --write-apidb host="localhost:5432" database="openstreetmap" user="osm" password="osm" populateCurrentTables=yes validateSchemaVersion=no

JOSM connects to local OSM server and download/upload data to local database. The problem is when we want to export Rails port database into OSM file to save updates and changes created by team members using the following line:
osmosis --read-apidb host="localhost:5432" database="openstreetmap" user="osm" password="osm" --write-xml file="from_pg.osm"

Executing the command throws the following error, how it should resolve?
org.openstreetmap.osmosis.core.OsmosisRuntimeException: Database version mismatch. The schema contains unexpected migrations [20170222134109, 20161011010929, 20150111192335, 20140115192822, 20191120140058, 20160822153055, 20140117185510, 20181020114000, 20180204153242, 20150110152606, 20190518115041, 20150818224516, 20190702193519, 20150222101847, 20190623093642, 20161002153425, 20131212124700, 20140519141742, 20140210003018, 20181031113522, 20190716173946, 20140507110937], may need to upgrade osmosis or specify validateSchemaVersion=no.
    at org.openstreetmap.osmosis.apidb.v0_6.impl.SchemaVersionValidator.validateDBVersion(SchemaVersionValidator.java:118)
    at org.openstreetmap.osmosis.apidb.v0_6.impl.SchemaVersionValidator.validateVersion(SchemaVersionValidator.java:53)
    at org.openstreetmap.osmosis.apidb.v0_6.ApidbReader.runImpl(ApidbReader.java:74)
    at org.openstreetmap.osmosis.apidb.v0_6.ApidbReader$1.doInTransactionWithoutResult(ApidbReader.java:102)
    at org.springframework.transaction.support.TransactionCallbackWithoutResult.doInTransaction(TransactionCallbackWithoutResult.java:34)
    at org.springframework.transaction.support.TransactionTemplate.execute(TransactionTemplate.java:133)
    at org.openstreetmap.osmosis.apidb.common.DatabaseContext2.executeWithinTransaction(DatabaseContext2.java:94)
    at org.openstreetmap.osmosis.apidb.v0_6.ApidbReader.run(ApidbReader.java:99)
    at java.lang.Thread.run(Thread.java:748)