texts: How to correctly install QGIS with plugins and GRASS-integration from source?

texts:
I'm running a Sabayon GNU/Linux (based on Gentoo) at my home desktop system and I want to install Quantum-GIS.

I'm currently using the qgis-1.7.0-package from the repositories. But this is somehow very minimal. It does not support downloading plugins (Plugins > Fetch Python Plugins disabled: I think this is to avoid security risks.) and it doesn't seem to be connected with GRASS in any way (at least in any visible way).

The official download guide is not very helpful as it only describes ways to install qgis in major linux distributions (Ubuntu, Debian, etc...). Anyways, this guide suggests to look out for packages like python-qgis or qgis-plugin-grass. This seems to be what I am looking for, but it is not included in Sabayon/Gentoo-repositories. (This seems to be a major issue with any non-Ubuntu/non-Debian Linux-distribution.)

My question is, how to install Quantum-GIS with full python-plugin-support und full GRASS-plugin-integration from source, where to get required source code for everything and how to compile it correctly?

The result should look like something I've found in this comment pointing to this video tutorial.

Update 01/01/2013: Question now focuses on compiling all packages on my own. I found out repositories are not very helpful with this issue [1,2,3,4,5].


texts: How to synchronize the marker's rotation field with the angular distortion in the map projection

texts:
I'm using QGIS 1.8.0 to plot a global dataset showing the positions of aircraft in flight.

The data I have are .csv files with three columns: Latitude, Longitude, and angle of travel.

I assign the "angle of travel" column of data to the rotation field in the "Style" tab of the "Layer Properties" dialog window.

This works fine except that that the rotation angle does not take into account the angular distortion of the map projection.

This means that in places where the angular distortion is high, the flight path looks like a row of dominoes:

Whereas where the angular distortion is low the flight paths looks like a straight line, which is what I want.

So my question is:

In QGIS is there a way to synchronise the marker's rotation field (now called "angle" in QGIS 2.0.1) with the angular distortion in the map projection?

EDIT:

Since posting this question I have started using the new QGIS 2.0.1.
In QGIS 1.8.0 I don't think it was possible to do what I am asking, but in the new version it looks like it is. However, it requires a thorough understanding of the maths behind map projections, which I do not have.

In the "Style" tab of the "Layer Properties" dialogue window, when assigning "data defined properties" to the angle of the marker, it is now possible to open an "Expression string builder" dialogue window.

Here, I presume, it will be possible to enter a mathematical equation which modifies the angle of the marker so that it accounts for the angular distortion in any particular map projection.

The map projections I am interested in are mollweide (or other similar elliptical projection) for a global overview,  and also the UTM zones.

The reasons for my map projection choices are:

Mollweide:
I like the look of mollweide for showing global trends.

UTM:
To make a complete visualization using all the data and revealing some of the lesser used routes clearly I will have to import and rasterize the data in geographic portions as I don't have enough RAM to deal with the entire dataset in one go.

I was thinking of using UTM zone strips, a bit like a set of globe gore, which I can then stitch back together once rasterized using Geocart.

The Mercator projection idea below is not ideal as it creates a pronounced narrowing of features near the poles, and there is a fair amount of air traffic in the north.
The conformal projection idea represents a lot more work for me, as many more tiles will be needed than if I use pole to pole UTM strips.

So, I'm still holding out for someone to tell me at least whether or not it's theoretically possible to type in a formula that will affect all the points in the correct way, if not the formulae themselves.


texts: Generating NDVI with multi channel images?

texts:
I have a 6 channel proprietary camera and I need to create a NDVI map. I created a map from each band, so now I have 6 files each one on a different band width.

How can I proceed now using QGIS?


texts: How to convert from lat,lng to tile coordinates?

texts:
Given a pair of lat/long coordinates (and a screen size if needed), we need to establish the Z/X/Y coordinates for the map tiles which will be requested when we tell our map fitBounds(coordinatePair);.

Background

We're trying to load test our map application. Our test partner uses NeoLoad, which I don't understand too well, but I think it stands in place of the page on the browser.  It seems to fire requests based on the recording, and allows us to parameterise them to simulate real usage - rather than repeatedly doing exactly the same thing (and getting the same piece of land in the resulting map).

We're having trouble parameterising the tile requests, which is unfortunate, because that's the main purpose of the load test.

Our map allows users to enter a text search string to help them find a piece of land.  Then we figure out the coordinatePair and show it on the map.  But while coordinatePair is in degrees, the resulting requests to the map server use the Z/X/Y coordinate system.

So we can't get our load tests to get map tiles corresponding to the right piece of land.


texts: User memory limit exceeded when applying iterative procedure to am image collection

texts:
Please help me understand earth engine memory usage. I am trying to apply an iterative procedure to a collection of images and then combine these images on a per-pixel basis. Here is the code structure that I am using:
var step = function(iteration, list) {

  var h0 = ee.Image(ee.List(list).get(-1));
  var h = h0.resample('bilinear').unmask(999);
  ...
  manipulations with h, involving several temporary variables being created
  ...
  return ee.List(list).replace(h0,h); // make sure not grow the list with intermediate steps
};

var iterate = function(h0) {
  var first = ee.List([h0]);
  var seq = ee.List.sequence(1,40);
  var tmp = ee.List(seq.iterate(step,first));
  return ee.Image(tmp.get(-1));
};

var G12 = imageCollection.filterMetadata('grid_id',"equals",12); //collection of uploaded images
var R12 = G12.map(iterate); // filtered collection: 99 images
var h1 = R12.toArray().arraySort().arrayPad([99,1]).arrayGet([90,0]);

So, function "step" defines one step of an iterative procedure. Having an old school coding background, I would assume that temporary variables created within the function are being pushed into the stack and then poped after the function has been executed (no memory is wasted). The list variable that is being passed in and out of the function does not grow, I merely replace the last element in the list. Function "iterate" creates a temporary list and maps the "step" function over it in iterative manner (each step using the result of the previous step). It takes an image as the input and spits out a processed (iterated) version of the image. The function was tested and works as intended.
Now I am trying to map it over an image collection. The result is being reduced to a single image via per-pixel sorting procedure. Now, I know, converting an image collection into an array is a memory hog, but I do not see any other way of doing this. More importantly, the last line calculation (reduction) works fine if the number of iterations in "iterate" function is below 40. Once it gets above 40 I get "user memory limit exceeded" error. This tells me that the iterative process adds to memory usage with each iteration, which I would not expect in a traditional code. Is my interpretation correct? Is there any way to code this differently?


texts: Creating a custom attribute table in qgis print composer

texts:
is it possible to arrange header and content of an inserted attribute table in qgis print composer? I want to create a custom datasheet with a small map and listed feature information as follows:

I thought to change header font into white (and insert a label for each attribute content) but apparently this is not possible in qgis 2.0.


texts:
What criteria do you need to take account to before you choose geometry or geography map?

Is there a book or page about it?


texts: Interactivity in CartoDb with createLayer method?

texts:
I have a problem with my interactivity and I can't solve alone. When mouse hover are on any point I would like show by console.log any column of my dataset.

With my code, only show me cartodb_id.

var map;
function init(){
  // initiate leaflet map
  map = new L.Map('cartodb-map', {
    center: [40,-98],
    zoom: 4
  })

  L.tileLayer('https://dnv9my2eseobd.cloudfront.net/v3/cartodb.map-4xtxp73f/{z}/{x}/{y}.png', {
    attribution: 'Mapbox Terms &amp; Feedback'
  }).addTo(map);

  var layerUrl = 'https://mcnb.cartodb.com/api/v2/viz/85aef586-5c57-11e5-9915-0e73ffd62169/viz.json';

  // change the query for the first layer
  var subLayerOptions = {
    sql: "SELECT * FROM mcnb_dev",
    cartocss: "#example_cartodbjs_1{marker-fill: #109DCD; marker-width: 5; marker-line-color: white; marker-line-width: 0;}"
  }

  cartodb.createLayer(map, layerUrl)
   .addTo(map)
   .on('done', function(layer) {

    layer.on('featureOver', function(e, pos, pixel, data) {
        // print data to console log

            console.log("Event #" + data.cartodb_id + ", Name " + data.name + ", Clase: " + data.class+ ",Incertidumbre(m.) " + data.coordinateuncertaintyinmeters );
          });

      // change the query for the first layer
      var subLayerOptions = {
          interactivity: "cartodb_id,class,coordinateuncertaintyinmeters",
          sql: "SELECT * FROM mcnb_dev",
          cartocss: "#example_cartodbjs_1{marker-fill: #109DCD; marker-width: 5; marker-line-color: white; marker-line-width: 0;}"
            }

      var sublayer = layer.getSubLayer(0);
      sublayer.set(subLayerOptions);

     // sublayer.infowindow.set('template', $('#infowindow_template').html());

      sublayer.on('featureClick', function(e, latlng, pos, data) {
              alert("Hey! You clicked " + data.cartodb_id);
        });

      }).on('error', function() {
            //log the error
    });
}

SubLayer featureClick works but featureOver not.


texts: Visualize classes from shapefile attribute table

texts:
I have a shapefile with classification and I would like to see that map in QGIS.
When I open it, adding it as layer, I just see the overall area filled with the shapes  of the various classes but I have not any legend available or a QGIS layer setting file. How can I set a legend?


texts:
I am trying to add simple legend to the map.

I wrote some simple test code below to see what happens but it did not work.

var map = activeView.FocusMap;

            IGraphicsContainer graphicsContainer = activeView.GraphicsContainer;
            IMapFrame mapFrame = graphicsContainer.FindFrame(map) as IMapFrame;
            UID uid = new UID();
            uid.Value = "esriCarto.Legend";
            IMapSurround mapSurround = map.CreateMapSurround(uid, null);
mapFrame.CreateSurroundFrame(uid, null);

            var legend = (ILegend)mapSurround;

            var legendItem = new VerticalLegendItemClass();
            legendItem.HeadingSymbol.Text = "served";
            legend.AddItem(legendItem);

            var legendItema = new VerticalLegendItemClass();
            legendItema.HeadingSymbol.Text = "unserved";
            legend.AddItem(legendItema);

            IElement element = legend as IElement;
            element.Geometry = envelope as IGeometry;
            graphicsContainer.AddElement(element, 0);

            activeView.Refresh();

The legend I am after, must be similar to this:

Can anyone assist?


texts: I'm fairly new to QGIS, and I'm trying to create a buffer around the outline of Spain (including the coastline and boarders). When I do so, it buffers all the country from center to the edge, I just need the edge if that makes sense.


texts: Why are SRTM tiles (*.tif) not visible in my QGIS project?

texts:
I am wondering, why my SRTM tiles are not displaying correctly in my QGIS (3.10.4 LTR) project. If I load vector layers to an empty QGIS project first, they are visible and the later added SRTM are not visible. Other raster layers are visible.

The problem exists also vice versa: If I load the SRTM tiles into an empty QGIS Project first, they are displayed correctly. But if I load vector layers later to that project, the vector layers don't show up...


texts:
I have been trying to create contours using the attached XYZ file.
I keep getting an error.
What could be wrong?
I have the XYZ file, I can't seem to attach it here
The coordinates are in "Projected Coordinates (Africa) Nigerian West Best Projection".


texts: GDAL Translate GeoTIFF to ASC creates a mirror image

texts:
I created a heat map using the HeatMap plugin in QGIS from spatial points with EPSG: 4326 CSR. I want to use this layer as a bias layer for MaxEnt which only takes files in ASC format. However, Heat Map doesn't have the option for creating density maps in ASC directly so my output is in GeoTIFF. I am trying to translate the GeoTIFF to ASC format but GDAL is creating a mirrored image of my layer. The top half of the picture is the original final.

I am not sure what I am doing wrong. Is there an alternative procedure to translate my GeoTIFF? Or is there a way to fix the mirrored layer?


texts: G_realloc: unable to allocate 46320000 bytes of memory while managing large .xyz through v.in.ascii

texts:
I tried to convert large (~180MB) .xyz file into vector through v.in.ascii in QGIS Brighton (2.6.0) Grass plugin (it should be Grass version 6.4.3). The "primitive registration" stage lasted very long and then ended with following error: G_realloc: unable to allocate 46320000 bytes of memory at cindex.c:108. It seems that GRASS can't allocate more than ~45MB of input data; 45MB is tiny compared to my RAM of 2GB, so the problem is probably caused by some limit in GRASS or in QGIS, and probably related specifically to v.in.ascii method, because r.in.xyz worked and didn't throw this problem.

Why is there such a limitation, and how to bypass it?

According to a comment in a related question, this problem should be fixed in Grass 7. However, I wish to use GRASS in the QGIS plugin form, and it will take months before it will be upgraded to GRASS 7.

Is there any way how to solve it with GRASS 6?

OS: Windows 7 (64-bit)
QGIS architecture: 32-bit
SWAP file: 4096MB

I doubled the SWAP file size and the error is the same.


texts: Iterating two different feature classes based on sequential order and clip

texts:
I have a watershed feature class with 273 polygons, in the attribute table there is a SiteID column that is populated as so:
A1
A2
A3
A4
etc.

I also have a buffer feature class with 273 polygons with a SiteID column populated as so:
B1
B2
B3
B4
etc.

I want to use model builder to iterate through each of the features classes and to clip Watershed polygon labeled A1 with Buffer polygon A1, and the Watershed B2 with Buffer B2 and so on.

I have tried model builder but couldn't get it to run properly, I tried batch processing after splitting both feature classes based on their SiteID column but it wont populated the input and clip parameters properly..

How can I make this happen?

I searched but did not find this exact question here on GIS Stack Exchange.

ArcPy 10.0 iterate through shapefile1 clip shapefile2 to shapefile1's polygons

How to Batch Clip in ArcGIS Desktop 10 using Python/ArcPy?


texts: Adding spatial index to spatialite file in python using OGR?

texts:
I have a spatialite data set without a spatial index on the geometry column and I would like to create it.
In OGR command line I can use the following to create it (OGR Spatialite driver):
ogrinfo db.sqlite -sql "SELECT CreateSpatialIndex('the_table','GEOMETRY')"

I would like to do this using the GDAL/OGR python bindings instead and I have tried to mimic the command line procedure:
driver = ogr.GetDriverByName('SQLite')
data_source = driver.CreateDataSource('db.sqlite', ['SPATIALITE=YES'])
layer = data_source.CreateLayer('the_table', None, ogr.wkbLineString25D, ['SPATIAL_INDEX=NO'])

data_source.ExecuteSQL('SELECT CreateSpatialIndex(\'the_table\', \'' + layer.GetGeometryColumn() + '\'')

The above returns the error message:
ERROR 1: In ExecuteSQL(): sqlite3_prepare_v2(SELECT CreateSpatialIndex('the_table', 'GEOMETRY'):
near "'GEOMETRY'": syntax error

What should I do different to make this work?
The answer to the obvious question, "why do I create a spatialite file without the spatial index just to create it afterwards?" is because I need to do bulk insertion.


texts:
I installed Qgis Master  64 bits on win7 64 16gb of ram and I am unable to locate the georeferencer. All other Gdal application are there, Gdal is active, but no georeferencer. I tried yesterday, tried also reinstall and also I anstalled again today. No georeferencer at all.
Any Ideas?


texts: How to create a polygon feature using Python?

texts:
I'm making a plugin for Quantum GIS in python. I have a polygon layer and want to create a new feature.

So far, I have this because I need a polygon geometry to add the feature:

surface = QgsGeometry.fromPolygon (polygon)

The problem is that I do not create a polygon. I have a list of points where I want to place the polygon.

Can anyone help?


texts: Using composite keys for table joins in QGIS3

texts:
How to define two fields for join tables in QGIS3?

Screenshot of the option layer properties >> Joins:

Today, in the layer properties >> Joins option, you can only add one field for reference in tables/layer.


texts: What data projection should I use for area analyses in different countries at a global scale

texts:
I'm creating A4 maps to show where forest cover and mining fields overlap across the globe. Each A4 map focuses on a specific country (e.g. India, USA, China, Australia...). For each country, I need to calculate how much overlap / intersection there is between the forest cover in that country and the mining fields to work out the area of forest that is threatened by mining activity.

I have sourced / received the following data:

World country boundaries (polygons, Geographic coordinate system WGS_1984, Datum D_WGS_1984)

Global forest cover (polygons, Geographic coordinate system WGS_1984, Datum D_WGS_1984)

Mining areas (polygons in multiple Geographic coordinate systems e.g.)
Australia - GCS GDA94
India - GCS WGS_1984
Indonesia - GCS WGS_1984

To calculate areas (including total mining area, total forest area and overlap between the two) I need to work in a projected coordinate system (as area calculations are disabled in global coordinate systems).

Can someone please advise which projection I should chose?

Normally I would opt to work in a UTM zone but when dealing with large countries like Australia this warps the shape of the country on the map.

The area calculations (in sqkm) are also altered quite dramatically depending on the projection I use.

Can someone recommend a projected coordinate system that I can suitably use for all the country maps and that will produce the most accurate area calculations?


texts: Finding centre point of a raster in ArcGIS?

texts:
If I have, say, a rectangular and georeferenced raster file loaded in an ArcMap 10 document (.tiff w/ associated tfw), how do I easily find its centre point and store that point in a points vector layer?

Also, if I have multiple such rasters in my ArcMap document, how do I apply the process to all of them?

Unfortunately, I have zero Python experience. Therefore a programmatic solution is OK, but I will need specific instructions on how to load an existing script into ArcGIS 10, and run it on the rasters in questions. (BTW, the rasters are all in their separate layers).


texts: Optimizing combining multiple feature classes into one via arcpy

texts:
I have about 1000 feature classes spread over about 50 geodatabases which I am trying to append into an empty layer in a separate geodatabase.

When I previously did this on shapefiles, the process completed quickly. However due to earlier processes (requiring Query Tables) I had to move the data into a GDB, and this attempt to append the layers has been running for about 50 minutes.

Are there any recommendations to decrease the time required to merge these feature classes?  Should I copy the feature classes to shapefiles first?

The layers all have the same field set; I am not using any field mappings.

General process:

Loop through the GDBs using arcpy.da.walk; add feature class path/name into a list
Create a feature class in the output GDB using one of the fcs above as a template
Use arcpy.Append_management(fclist, out_layername)

Note: While the first step (obtaining the list of FCs via arcpy.da.walk) takes a considerable amount of time, that is acceptable - I'd like to improve the append or merge step.

The specific code snippet I am trying to improve is:

# Set the workspace environment setting
arcpy.env.workspace = os.path.join(rootworkingdir, gdbName)
out_path = gdbpath
out_name = fcname
geometry_type = "POLYLINE"
template = fclist[0]
# Execute CreateFeatureclass
arcpy.CreateFeatureclass_management(
    out_path, out_name, geometry_type, template)
out_name = out_path + "/" + out_name
arcpy.Append_management(fclist, out_layername)

The total number of expected features is about 5000.


texts: Randomly generate points using weights from raster

texts:
I need to randomly generate ~ 40K points throughout my region of interest, but I need to weight the point selection by a set of probability weights, stored as a raster  - such that higher probability areas have more of the random points and lower probability areas have fewer points. Additionally, I need to indicate a minimum distance between the points.
I found that I could do this in ArcMap using the "create spatially balanced points" tool in the GeoStatistical Analyst package - but I do not have a license for that package.
Is there a way to do this in QGIS or R?
 I found the sp package in R to randomly select points - but I do not see how I can weight that selection given the "type" options.

Any suggestions on how to solve this problem? Other R packages that I should check out?


texts: Creating layer for text labels in QGIS 3.10

texts:
How can I create a layer to hold individually placed text labels in 3.10?
I can do it in print composer, but would like to have it as a layer for repeated use.
The labels are simple text such as  "Barecrofts", labelling locations and to be placed manually.


texts: Map To KML does not work with mxds that contain service layer

texts:
I am trying to use the map to kml conversion tool in a python script. The mxd I provide as input, contains data coming from map service layers. I have no security enabled to that particular layers.
I also try to run the specific tool in ArcMap providing mxds that contain layers coming from a map service but I get the same error:

Unable to open MXD - The document contains secure (password-protected) layers.

If I run the same tool with an mxd that contains data coming from a local geodatabase, but not from a published service, it works properly.

Can anybody tell me why the map to kml tool does not work with map service layers?


texts: Finding the LineString Corresponding to a Certain Road in a Roads Database

texts:
I am using tiger line road data to find roads that runs through certain points. I know what road is immediately next to the point by the road name and road type. What I currently do is grab all the road segments within a 10km radius from the point of interest and then filter those based on the name of the initial road segment and its type. However, road names change and therefore, I fail to always create the correct road that I need.
SELECT [fields] FROM point_and_initial_road_matches pr INNER JOIN tl_edges tl ON ST_DWithin(ST_GeomFromText(pr.point_wkt, 4326), tl.geom, 0.1) AND SIMILARITY(pr.road_name, tl.fullname) > 0.7 AND pr.road_type=tl.mtfcc

If you look at the image below, although I am trying to catch the highway, because the name of the highway changes, I only am able to extract part of the highway and the rest is coming from a road that is next to it which has a similar name. How can I find line strings that are connected, so that, I can restrict the set by the road type and extract only the line strings that have the same road type as the original line segment? or is there a better way to do this whole exercise?


texts:
I am trying to reclassify a DEM by using  the new Reclassify by layer tool in QGIS3. I am using the following input table:

However, as shown in below screenshot of the tool, the Minimum class value field and Maximum class value field don't show any values.

What could I have done wrong here?


texts: Combining DEM raster with vector data in ArcScene?

texts: Ive have a raster layer of earthquake data that includes depth.  I also have a DEM file.  The problem is that I need to plot earthquake depth relative to the DEM file rather than sealevel in order to realistically place the focal points in ArcScene 9.3.   Simply using the depth data for base heights results in earthquakes that occur in the water rather than below the sea floor.  I've tried fooling with the Z unit conversion and offset, but any adjust that produces reasonable result for earthquakes beneath the ocean floor leaves earthquakes that occur under mountainous area looking too deep.


texts: mbtiles-extracts does not work with my geojson

texts:
I am working with the mbtiles-extracts library: https://github.com/mapbox/mbtiles-extracts. Given an mbtiles file and a geojson, the library extracts mbtiles cutouts of the polygons in the geojson.

Except, it does not work for some geojsons.

As documented in https://github.com/mapbox/mbtiles-extracts/issues/10, it does not work for any polygon present in the following geojson for the mbtiles listed below:

MBTILES: Pretty much every tileset on https://osmlab.github.io/osm-qa-tiles/historic.html. Can definitely replicate it using for https://s3.amazonaws.com/mapbox/osm-qa-tiles/historic/2008-Q1-qa-tiles.mbtiles
GEOJSON: https://gist.github.com/arunasank/c98ef28abd51b54579c7609f8ebc1de1

When I debug, I see that:

This function: https://github.com/mapbox/mbtiles-extracts/blob/master/index.js#L175
returns a falsy value for all the features in the above geojson, resulting in a falsy value here: https://github.com/mapbox/mbtiles-extracts/blob/master/index.js#L58-L62

As suggested by Nikhil V J on https://lists.openstreetmap.org/pipermail/talk-in/2018-November/003211.html, I also tried validating the geojsons, and ensuring that they follow the right hand verification rule, but that did not fix things.

Any help in understanding why these polygons are different from other polygons, or if there's a pattern to these polygons that I am missing, would be really helpful!


texts: OpenStreetMap tiles not loaded despite being received

texts:
While using OpenLayers to load OpenStreetMap, I observe that part of the map would sometimes be covered with blank pink tiles. I searched and found similar complaints with solutions suggesting the modification of OpenLayers.IMAGE_RELOAD_ATTEMPTS parameter to an integer number greater than 0. I did exactly that, but it doesn't solve the problem.

Using Firebug, I am able to verify that the tiles were indeed downloaded successfully. However, for some unknown reason, they just do not load onto the map. I have no issue with Google Maps as a base layer. The problem seems to be restricted to OpenStreetMap only. Any idea how to solve this pink tile of death?


texts:
According with this link
http://postgis.net/workshops/postgis-intro/advanced_geometry_construction.html
I'm creating lines that describe routes.
I have a table stopcentroids4 like in figure with the following fields:
gid, routes, geom

and i want to draw a line for each route.

Adjusting the query with my table name (see the code beloww) as described at the link above i get that the query runs but with a wrong output.
For example the line drawn for routes='1' connects all stopcentroids points whose associated routes contain the character '1' (that are also 1A, 16, 21B and so on) as you can see in the picture with the wrong route..

I want the route that connects points associated only to the same character '1' (and not points associated to 11, 21, 12A and so on)as shown in the picture with the right route.

How do I correct this query?

This is my function

> CREATE OR REPLACE FUNCTION walk_centers(bigint, text)   RETURNS
> geometry AS $BODY$ WITH RECURSIVE next_stop(geom, idlist) AS (
>                  (SELECT geom AS geom,
>                     ARRAY[gid] AS idlist
>                    FROM stopcentroids4
>                   WHERE gid = $1)
>         UNION ALL
>                 ( SELECT s.geom AS geom,
>                     array_append(n.idlist, s.gid) AS idlist
>                    FROM stopcentroids4 s,
>                     next_stop n
>                   WHERE strpos(s.routes::text,$2) != 0 AND NOT n.idlist @> ARRAY[s.gid]
>                   ORDER BY st_distance(n.geom, s.geom) ASC
>                  LIMIT 1)
>         )  SELECT ST_MakeLine(geom) AS geom FROM next_stop; $BODY$   LANGUAGE sql VOLATILE   COST 100; ALTER FUNCTION walk_centers(bigint,
> text)   OWNER TO postgres;
>
> This is my code
>
>
> CREATE TABLE lines AS
> -- Distinct route identifiers! WITH routes AS (   SELECT DISTINCT unnest(string_to_array(routes,',')) AS route   FROM stopcentroids4
> ORDER BY route ),
> -- Joined back to stops! Every route has all its stops! stops AS (   SELECT s.gid, s.geom, r.route   FROM routes r   JOIN stopcentroids4 s
> ON (strpos(s.routes, r.route) <> 0) ),
> -- Collects stops by routes and calculate centroid! centers AS (   SELECT ST_Centroid(ST_Collect(geom)) AS geom, route   FROM stops
> GROUP BY route ),
> -- Calculate stop/center distance for each stop in each route. stops_distance AS (   SELECT s.*, ST_Distance(s.geom, c.geom) AS
> distance   FROM stops s JOIN centers c   ON (s.route = c.route)
> ORDER BY route, distance DESC ),
> -- Filter out just the furthest stop/center pairs. CREATE TABLE brighton_lines AS
> -- Distinct route identifiers! WITH routes AS (   SELECT DISTINCT unnest(string_to_array(routes,',')) AS route   FROM stopcentroids4
> ORDER BY route ),
> -- Joined back to stops! Every route has all its stops! stops AS (   SELECT s.gid, s.geom, r.route   FROM routes r   JOIN stopcentroids4 s
> ON (strpos(s.routes, r.route) <> 0) ),
> -- Collects stops by routes and calculate centroid! centers AS (   SELECT ST_Centroid(ST_Collect(geom)) AS geom, route   FROM stops
> GROUP BY route ),
> -- Calculate stop/center distance for each stop in each route. stops_distance AS (   SELECT s.*, ST_Distance(s.geom, c.geom) AS
> distance   FROM stops s JOIN centers c   ON (s.route = c.route)
> ORDER BY route, distance DESC ),
> -- Filter out just the furthest stop/center pairs. first_stops AS (   SELECT DISTINCT ON (route) stops_distance.*   FROM stops_distance )
> -- Pass the route/stop information into the linear route generation function! SELECT   row_number() over () AS gid, -- QGIS likes numeric
> primary keys   route,   walk_centers(gid, route) AS geom FROM
> first_stops;
>
> -- Do some housekeeping too ALTER TABLE lines ADD PRIMARY KEY (gid);


texts: Creating Mosaic Dataset in ArcMap with Colormap as renderer or adding Colormap after creation?

texts:
I have several rasters in tiff format that have a colormap assigned to them:

I then create a Mosaic Dataset with specs matching the imagery and add all of my tiffs to it (you can ignore the warnings this is just to illustrate what i did):

Colormap does not show up and it does not let me import my Colormap from the tiffs:

Is there any way to add/import my Colormaps so that we can view the data according to those values?


texts: Store live tracked GPS points with additional attributes (alitude, DOP, accuracy)

texts:
I am using the Live GPS tracking feature from QGIS (GPS Information Panel).

For a project I need to collect points in the field and directly see and store them in QGIS. For this I have a setup of a GNSS-Receiver (Trimble R1) and a Windows 10 tablet that I can connect with Bluetooth. On the tablet I installed QGIS 2.18.2 Las Palmas. The receiver is able to connect to the QGIS GPS-Information Panel. It is also possible to “add points”. However, if a point is added, it is no longer possible to access the additional GPS information (like DOP or accuracy).

Is it possible to save points using live GPS tracking that include all the GPS Information (time, altitude, different dop, accuracy) as attributes, and if yes, then how can it be done?

I do not think that this is a duplicate of Add "Position and additional attributes" in table because the answer there does not give me enough clarification.

UPDATE:

I decided to use following bumpy workaround to tackle this problem for now:

First I make sure that the GPS information panel writes a NMEA log file while being in the field. Points are stored at the required positions using the normal “add point” feature tool of GPS information panel. Back home NMEA2QGIS2 plugin is used to load the NMEA log file into QGIS. The GPS points are buffered and then spatially intersected with the closest record point from the NMEA log file.

It might not be the most precise way to do this, at least it helps me to save details about the GPS points, which help to judge the quality of the GPS also in retrospective. If somebody comes up with a better solution, I’d be very happy to hear about it!


texts: Which Geodatabase supports SQL and is easily portable?

texts:
I want to choose the geodatabase type, which could be:

permits direct sql queries.
easily portable, because the app would installed in one computer.

I know that the personal geodatabase is a good choice, but I would like to know if there's some limitations in personal geodatabase apart the size (the size won't be a problem).

The geodatabase would be used for some geoprocessing operations and some sql queries.


texts: Can I use a Geotools raster layer in a Geomajas webapplication?

texts:
I am working on a Geomajas web application and have a Geotools RasterLayer (a GridReaderLayer to be precise) that needs to be integrated with the map in the Geomajas app.

I have been searching the web for quite a while now but didn't find anything really useful. Can this even be done? And if yes, how do I get started?

I would be extremely grateful for any pointer in the right direction!

Wouter


texts:
I am new to GEE. I have an application where values can be set then points digitized.  I need to programmatically apply the values to fields(properties) while a user is adding points. I don't seem to find an API into the 'Configure geometry import' tool where the attribute field values can be set then applied to any new points that are digitized. I don't want a feature collection because each feature can have different values for the limited set of fields.

Anyone have an example of digitizing points and applying different values at different times while to user is adding points?


texts: Inserting C-ordinates (LAT. and LONG.) in map composer of QGIS

texts:
In QGIS 2.8, how do I insert co-ordinates (degree, minute and second) in map composer?

When I did it, the the X and Y intervals value were showing on the frame instead of actual lat. & long of the area. Intervals unit were set as Map Units. co-ordinate Format as Degree, Minute and second with suffix. But the coordinate appearing on the canvas is not correct.


texts: OpenLayers/Geoserver rendering transformation is failing at the Prime Meridian

texts:
I have a layer that represents surface pressure in GeoServer. I am using a rendering transformation SLD to generate the contours. When I try to map it in OpenLayers, the contours show up fine when zoomed out all of the way.

However, when you zoom in, the contours disappear at the Prime Meridian. I'm thinking there may be an issue with how GeoServer or the SLD are handling this area and am wondering if anyone knows of a solution.

I'm not totally sure if this is an issue within OpenLayers or if it is an issue with my SLD in GeoServer. I'm thinking it is in GeoServer specifically with the rendering transformation because layers that do not have a rendering transformation show up fine.

Please see a screenshot of the issue below.


texts: Using population data as Z-axis to make 3D map in Python?

texts:
I want to make a 3D thematic map with the different status' population data as the height in Python. I have the shapefile of the status and the population data, although I don't know how to set the population attribute data as the Z-axis. I want to make a map like the following.

Is there a python package that allows the 3-D mapping of shapefiles such that I can specify an attribute as the Z value?


texts: Merge thousands of rasters with gdal_merge.py

texts:
I have several thousands of small rasters, kb size, tens or hundreds of pixels most of them. Batch command of gdal_merge.py can't take all the files at once, so I am doing partial merges of chunks of 500 images, and that's relatively fast.

But when I try to merge the partial merges it takes way too much. In one even gdal_merge.py gave up and returned Killed.

Has anyone experienced something similar, and found an optimization?

Alternatives to gdal_merge.py are welcome


texts: Debugging ERROR 000967 from standalone script?

texts:
The following script is working in the Python window of ArcMap but not standalone in the PythonWin IDE:

   import arcpy
   from arcpy import env

   env.workspace = r"J:\Non_RS_MXDs\LongHaul.gdb"

   mxd = arcpy.mapping.MapDocument(r"J:\Non_RS_MXDs\RSF_Long_Haul.mxd")
   df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
   sourceLayer = r"J:\Non_RS_MXDs\Layer Files\Facility.lyr"
   for i in df:
       arcpy.ApplySymbologyFromLayer_management(i, sourceLayer)
   print "Done!"

When I run it I get this error:

    ERROR 000967: The input layer type must be a feature layer, raster layer, or tin layer.
    Failed to execute (ApplySymbologyFromLayer).


texts: Warping a Polar Stereographic image to UKNG projection

texts:
I have a North Polar Stereographic image over the UK. The top left and bottom right coordinates are: (-17.644 E, 59.683 N) (9.231 E, 45.25 N).
I want to reproject this to UKNG (aka BNG) projection.
This is my code. First I using gdal_translate to append georeferencing data to my png, output as tif.
PROJ="+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs"

gdal_translate -of GTiff -a_ullr -17.644 59.683 9.231 45.25 -a_srs "$PROJ" the_image.png test.tif

Alright, so now I have geotiff with the appropriate bounds, I want to warp it to the UKNG projection
UKNGPROJ="+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs"

gdalwarp -s_srs "$PROJ" -t_srs "$UKNGPROJ" test.tif plswork.tif

I then convert back to png so I can overlay on my OpenLayers map:
convert plswork.tif ohwhydoesntthiswork.png

And I get dodgy output. The image is diagonal and doesn't overlay onto the map correctly at all.
Both NP Stere and UKNG proj definitons come from the EPSG site, https://epsg.io/27700 , https://epsg.io/3411


texts: What is the best way to manage dependencies between layers?

texts:
I have multiple layers that has dependency with each other.
Such as: City; Courts; Lots; Front Lot; Buildings.
The polygon's "block" depends on the polygon's "City." Example dependency: Buildings-> Front Lot -> Lots -> Blocks -> City.

On my server I have installed the Oracle SDE.
He has the ESRI SDK. With the "spatial functions" for Oracle.

I thought I'd solve this with the "spatial functions."
I think to check the intersection of the layers. As the example of the link:

http://resources.arcgis.com/en/help/main/10.1/index.html#/A_quick_tour_of_SQL_functions_used_with_ST_Geometry/006z0000003n000000/

In my future solution, when a building is in another district must update this value automatically.
Maybe even cut off part of a neighborhood that is outside the city limits.

What better way to manage these dependencies layers?
There is a ready tool? Or API? ADF/Python etc? Or another solution ready?

Thank you!


texts: When to store geometries using SQL-server geometry datatypes or ArcSDE?

texts:
With the release of SQL Server 2008 Microsoft added geometry data types. So geometries can now in an ESRI enviroment be either be stored using SQL-Server datatypes or ArcSDE datatypes.

A colleague of mine asked today when should we choose one product over the other?

We already have business databases in SQL-Server and in ArcSDE, if we need to create a new database, what options do we have? Are there any advantages/disadvantages to storing data in one way or other?

Refrased, Old Question is here:

"With the release of SQL Server 2008 Microsoft added geometry data types. So geometries can now either be stored in SQL-Server or ArcSDE.
A colleague of mine asked today when should we choose one product over the other?
We already have business databases in SQL-Server and in ArcSDE, if we need to store a new polyline, is there any reason to promote one place over the other? "


texts: How to control the arrow length of vector field markers in the layer list legend?

texts:
In QGIS 2.14.5 LTR (and recently in QGIS 3.0) I've set up a beautiful map design with a graduated renderer and vector field markers for each value range (arrow length is directly proportional to stream velocity, which ranges from 0.001 -~5m/s).

This worked great, but I wonder if I missed a setting to shorten the arrows in the legend. With a 'normal' width layer window, the class values cannot be seen:

Units: millimeter

My first guess was, that it has sth. to do with switching arrow length, line width units from millimeter to map units, but this doesn't affect the legend:

Units: map units

UPDATE

Issue still persists in QGIS 3.0.

I created a screencast to visualize this: https://youtu.be/oxJ3tGXFB18

In properties / symbology the symbols render correctly:

Note that the layer list displays the markers correctly, when duplicating the layer or applying changes to the style as long as layers panel width does not change.


texts: openlayers 4 method forEachLayerAtPixel() throws SecurityError: The operation is insecure

texts:
Working with GeoExt 3.1.0, ExtJS 6.2.0 and Open Layers 4.4.1.

I am trying to select a layer on my map with the singleclick event so I can use layer.getSource().getGetFeatureInfoUrl and display the information contained in the table row corresponding to the selected/clicked geometry on the layer (just like the layer preview option given by GeoServer in OpenLayers format). So far:

                map.on('singleclick', function(evt) {
                    document.getElementById('nodelist').innerHTML = "Loading... please wait...";
                    var view = map.getView();
                    var viewResolution = view.getResolution();

                    var pixel = evt.pixel;
                    //pixel = map.getEventPixel(evt.originalEvent); <-- tried with and without this line
                    map.forEachLayerAtPixel(evt.pixel, function(layer) {
                        var source = layer.getSource();
                        console.log(layer);
                        return true;
                    });
                    var url = source.getGetFeatureInfoUrl(
                            evt.coordinate, viewResolution, view.getProjection(),
                            {'INFO_FORMAT': 'text/html', 'FEATURE_COUNT': 50});
                    if (url) {
                        document.getElementById('nodelist').innerHTML = '<iframe seamless src="' + url + '"></iframe>';
                    }
                });

Problem:

The method forEachLayerAtPixel is causing this error SecurityError: The operation is insecure. to appear in the firebug console every time I click any point on the map (on or off any geometry)

I have read articles relating this problem to the Access-Control-Allow-Origin tag being missing, AFAIK this is resolved (in my case in GeoServer Jetty) by adding this chunk of code:

    <filter>
        <filter-name>cross-origin</filter-name>
        <filter-class>org.eclipse.jetty.servlets.CrossOriginFilter</filter-class>
    </filter>
    <filter-mapping>
        <filter-name>cross-origin</filter-name>
        <url-pattern>/*</url-pattern>
    </filter-mapping>

in:  C:\Program Files (x86)\Boundless\OpenGeo\jetty\etc\webdefault.xml

and including this {mode: 'cors'} along with the url being targeted in the AJAX call. At least this is how it is working for me.

Or maybe this has nothing to do with the problem, in which case this question just got a lot longer then it should have.


texts:
On this occasion I need someone to guide me, to implement cascading form controls in QGis 2.18. I am using a code for the attributes "area" "county" and "town" to perform a dynamic filter. The problem is, it does not filter the attributes. This is the code, I got it from https://medspx.fr/blog/Qgis/cascade_forms/ by Médéric RIBREUX.

Is it possible to change this code to work with QGIS 3?

# -*- encoding: utf-8 -*-

from PyQt4.QtCore import *
from PyQt4.QtGui import *

myDialog = None

def formOpen(dialog,layerid,featureid):
global myDialog
myDialog = dialog

# We need to introspect the form
# In QGis v2.4 autobuild form controls do not have a name... (fixed in v2.6)
# So we need to retrieve form controls by field order !
lst_children = dialog.findChildren(QComboBox)
area = lst_children[0]
county = lst_children[1]
town = lst_children[2]

# Clear the children QComboBox widgets
county.clear()
town.clear()

# When you change the value for area widget, change the content of the county widget
area.currentIndexChanged.connect(lambda: cascadeManage('COUNTY', 'ID_COUNTY', 'COUNTY', 'ID_AREA', county, area))
# When you change the value for county widget, change the content of the town widget
county.currentIndexChanged.connect(lambda: cascadeManage('TOWN', 'ID_TOWN', 'TOWN', 'ID_COUNTY', town, county))

def cascadeManage(layername, id_idx, txt_idx, fltr, widget, parent_widget):
'''Generic function to manage cascading QComboBox.
   * layername: the layer (flat table) to request values for this control
   * id_idx: primary key for the table (used in for the form control key)
   * txt_idx: field shown in the form control QComboBox
   * fltr: field which is used to filter the content of this widget based on the parent widget value
   * widget: widget object (QComboBox) to control
   * parent_widget: get the value from the parent QComboBox widget'''
from qgis.core import QgsMapLayerRegistry, QgsExpression

# Get the layer (flat table)
layers = QgsMapLayerRegistry.instance().mapLayersByName(layername)
if (len(layers) > 0):
    layer = layers[0]
else:
    return False

# get attributes indexes
ididx = layer.fieldNameIndex(id_idx)
txtidx = layer.fieldNameIndex(txt_idx)

# Get the value of the parent QComboBox
filter_code = str(parent_widget.itemData(parent_widget.currentIndex()))

# Build an expression with it
exp = QgsExpression("{0} = {1}".format(fltr, filter_code))
exp.prepare(layer.pendingFields())

# clear the current widget content
widget.clear()

# And fill it with filtered data:
for feature in layer.getFeatures():
    value = exp.evaluate(feature)
    if exp.hasEvalError():
        raise ValueError(exp.evalErrorString())
    if bool(value):
        widget.addItem(feature.attributes()[txtidx], feature.attributes()[ididx])


texts: Draw line between all points on a layer based on parent field

texts:
I have a table that contains georeferenced points, a bunch of other information and a parent for each point.  So let's say the unique identifier is name in the table, each record also has a parent name.

Using QGIS I want to create a layer containing a line between each point and its parent and be able to change the line colour based on some of the other information in the record.  I think I can do the line colour easily but automatically generating the table with the lines is something I am hoping someone can help me with.


texts: Enabling/disabling parameters based on Conditions/Preconditions in ModelBuilder?

texts:
I have been working with ModelBuilder and I am having trouble placing a precondition on an input parameter. However, I may be going about it in the wrong way and may just need a Python script.

My example will be based off of the 'Strip Map Index Features' tool (Cartography Tools --> Data Driven Pages --> SMIF) in ArcMap. When running the tool, the GUI shows a checkbox (Use Page Unit and Scale). There is also another parameter/textbox (Map Scale). The 'Map Scale' text box only becomes active or enabled when the checkbox is checked.

I am trying to recreate this functionality using ModelBuilder and I am not sure how to do so.


texts:
I am new to GIS and I wish to perform view shed analysis of wind turbines from a DEM.

I have downloaded OS Terrain 50 and the data is contained within numerous zipped folders.

I am using a Mac with OS X Mavericks and I would like to know how to unzip all the data and merge them into one file for QGIS?

Thanks in advance for any help.


texts: Python to read text file to add new field

texts:
So far I have written the following and both work:

arcpy.AddField_management("Wards", "test", "SHORT")

iFile =r'\file\Book1.txt'
f = open(iFile, 'r')

What I want to do is populate the new field 'test' with a field (value) from the iFile based on a match on two codes in the seperate file.

The text file has:

Code, Value

I think it might need updateCursor with a where statement ... or an if?

if (code in iFile == code in Wards):
    raise iFile(value)


texts: Interpolating surface as though bare earth DEM after using Extract by Mask in ArcMap

texts:
I have some stone wall features that I need to determine the height of in ArcGIS Desktop. I have extracted the features out of the DEM by using the extract by mask tool. Now I need a method for interpolating the surface where the features were removed so that the surface reflects bare earth with the features removed.

I tried to convert the clipped raster to points and then back to a raster using IDW interpolation but the process took a lot of time and ultimately failed.

Is there an easier method for doing this or am I maybe running into some issues with my approach?

Basically, I need to find the height of these stone walls by creating a bare earth DEM without them in it. I've tried to create a DTM from LiDAR data but the stone walls weren't completely removed. So after digitising them with polygons, I've been able to extract them out of the DEM and need to interpolate the surface to effectively create a DTM without the walls in so that I can accurately calculate their height.


texts:
I wish to select records based on the NRN attribute in my dataset :  "NRN"  IN ('1','2','23','31')

Some records have multiple values in the NRN field separated by the underscore character.
For example : 20_21, 14_23 etc

How do I select any records that contain the values that I want ignoring the records that have part of the string but not the whole number? Example : select 14_23 but not 20_21.


texts: How to add all tables from selected Personal geodatabase

texts:
I have a personal geodatabase which contains multiple tables and feature classes which are within feature datasets. I can add all the feature classes to ArcMap but I am unable to add all tables. Can anyone help me?

I just want to add all tables from personalGDB to ArcMap in single click


texts: How do I gather X and Y co-ordinate data from a vector grid overlay in QGIS 1.8.0?

texts: I have overlaid a 10m x 10m grid onto a UK Ordnance Survey raster map. I have set both raster and vector maps to the same CRS, (OSGB 1936 British National Grid; ESPG:27700). I have a csv file containing biological survey information gathered at differing 10m quadrats. I want to join the survey data to my 10m x 10m grid as a new layer. However, when I open the attribute table, only Y co-ordinate data is shown. How do I make QGIS show both x and y co-ord attributes for a grid so that I can import my tabular information?


texts: ESRI JavaScript API 3.8 Locator: SyntaxError {log: undefined, stack: (...), message: "Unexpected token <"}

texts:
I am having issues trying to set up a geocoder/locator on my server.  I keep getting an error SyntaxError {log: undefined, stack: (...), message: "Unexpected token <"}, but I can't seem to figure out why.  I tried exploring the error message, but I don't get much good information out of it.

I have tried using ESRI's geocoder widget then created my own set up using the locator service.  I have checked my address string, and I believe it should be passing correctly over to the server.

If I am understanding the error message right, it is saying there is a token "<" somewhere that shouldn't be there, but I cannot find it.

Other info...  I am running version 9.3.1 for my server, using the JavaScript API 3.8 that ESRI hosts, and created an address locator on my server using street centerlines for the area I am mapping.

Has anyone run into this issue before or have suggestions for me?  Right now, I am just trying to center the map on the location.  Once I get that done, I can worry about all the other stuff.

My code for the locator is below.Here is my development site.

    /////////////////////////////////////////////////////////////////
locator = new Locator(singleFieldGeocodeURL);
//add handler for completion
locator.on("address-to-locations-complete", showResults);

//trying to find the error
locator.on("error", function (error) {
    errorObject = error;
    console.log("error" + error);

});

//add handler for button click
on(dom.byId('locate'), 'click', function () {

    console.log(dom.byId('address').value);

    map.graphics.clear();
    var address = {
        "SingleLine": dom.byId('address').value
    };

    locator.outSpatialReference = map.spatialReference;

    var options = {
        address: { street: address },
        outFields: ["Loc_name"]

    };
    locator.addressToLocations(options);

});

function showResults(evt) {
    console.log("in show results");
    var candidate;
    var symbol = new SimpleMarkerSymbol();
    symbol.setStyle(SimpleMarkerSymbol.STYLE_SQUARE);
    symbol.setColor(new Color([153, 0, 51, 0.75]));

    var geom;
    arrayUtils.event(evt.addresses, function (candidate) {
        console.log(candidate.score);
        if (candidate.score > 80) {
            console.log(candidate.location);
            var attributes = {
                address: candidate.address,
                score: candidate.score,
                locatorName: candidate.attributes.Loc_name
            };
            geom = candidate.location;
            var graphic = new Graphic(geom, symbol, attributes, null);
            map.graphics.add(graphic);

            return false;

        }
    });

    if (geom !== undefined) {
        map.centerAndZoom(geom, 12);
    }

}


texts:
Attempting a run-of-the-mill ogr2ogr process with a KMZ file failed:

FAILURE:
Unable to open datasource `Rail.kmz' with the following drivers.

Is there a trick? (Obvious workaround is to unzip the .kmz first, but that's sometimes clunky, and you don't know what the internal structure will be exactly.)

GDAL 1.11.0, released 2014/04/16


texts: G_malloc error while running r.watershed in GRASS 7.0.4

texts:
I tried to use r.watershed function in GRASS GIS, and got an error message ERROR: G_malloc: unable to allocate 3837384596 bytes of memory at raster/r.watershed/ram/init_vars.c:146 before getting the outputs. Because I am not sure what caused this, so I would describe my situation from the beginning.

GRASS GIS 7.0.4 was installed by installing OSGeo4W package. (under Win 10, 64-bit, x64, RAM = 8 GB)
I created location 'nl' using GRASS GUI

set the EPSG code based on the actual location of the map I would load to GRASS GIS later

Then I set the extent and resolution of the location. (I loaded the map (later to be loaded in GRASS) to QGIS, and used the 'save as...' to check the extent of this map, and set the extent based on the information from QGIS.) But I didn't change the setting of resolution, just left them as default.

I created the mapset 'nl'

In 'nl' mapset, I used r.import to load the raster map, although there were some warnings, the result looked fine.

r.import input=C:\Users\Heinz\Desktop\GRASStest\hc.tif output=hc
WARNING: Datum <unknown> not recognised by GRASS and no parameters found
WARNING: Datum <unknown> not recognised by GRASS and no parameters found
WARNING: Datum <unknown> not recognised by GRASS and no parameters found
Proceeding with import of 1 raster bands...
Importing raster map <hc>...

I used r.watershed, but the error message occurred:

  r.watershed elevation=hc@nl threshold=5000 accumulation=acc
  drainage=dra basin=bas stream=str

  SECTION 1a (of 5): Initiating Memory.

  Current region rows: 57567, cols: 53957

  ERROR: G_malloc: unable to allocate 3837384596 bytes of memory at
  raster/r.watershed/ram/init_vars.c:146

  WARNING: Subprocess failed with exit code 1

  WARNING: category information for [acc] in [nl] missing or invalid

  WARNING: category information for [dra] in [nl] missing or invalid

  WARNING: category information for [bas] in [nl] missing or invalid

  WARNING: category information for [str] in [nl] missing or invalid

I don't know what stage-setting causing this error, and how to solve this error.

Update#1:
 I installed GRASS GIS 7.0.5 through OSGeo4W 64bit installer (under Win 10, 64-bit, x64, RAM = 8 GB), and ran r.watershed with the same raster file loaded into the mapset, but got the same error message I posted above.

It's weird because I successfully ran r.watershed on the same notebook before, just with GRASS GIS version 6. (This NB had been reformatted once.)

The result of g.region -p:

Besides, as I launched GRASS GIS GUI, some error messages showed up:

I am wondering if these errors in relation to the error message showed as I ran r.watershed, and I have checked that under C:\OSGeo4W64\apps\Python27\Lib\site-packages I got numpy installed.


texts: PostGIS raster MapAlgebra callback optimization

texts:
I'm trying to speed up an ST_MapAlgebra callback function that essentially takes a raster with four bands, checks each band for values less than a given number, and where the check is true, return the fourth band.

This is my messy function:

create or replace function
callback_fn(pixel float[], pos integer[], variadic userargs text[])
  returns float
  language plpgsql
  immutable
as $$
  declare
    roaddist float;
    protected_area float;
    citydist float;
    return_val float;
  begin

    citydist := pixel[1][1][1];
    roaddist := pixel[2][1][1];
    protected_area := pixel[3][1][1];
    return_val := pixel[4][1][1];

    if protected_area = 1 then return null;
    end if;
    if citydist <= 2000 then return return_val;
    elsif roaddist <= 500 then return return_val;
    else return null;
    end if;
  end;
$$;

So from what I understand, the function above will go through every pixel and either return the "return_val" or it will return NULL. So what I end up with is basically a single band raster with excluded areas removed.

This is exactly what I want, but I was hoping to run this faster. Is there a way to run multiple pixels simultaneously instead of going through each one individually? I'm confused why ST_Reclass can run so quickly to reclassify a whole raster while ST_MapAlgebra seems to take so much longer when (in my mind, at least) the functions are doing similar things - checking the values of a raster and where true, change value, otherwise set null.  Is there a way to use ST_Reclass in a callback function to speed up this process?

Additionally, is there a good source for documentation on writing PostGIS raster callback functions? I have been googling a lot over the last few days and have found extremely meager examples.


texts: Creating polygon feature class from X,Y coordinates using ArcPy?

texts:
I am trying to write a script that creates a new polygon feature class containing a single ( square ) polygon with the following coordinates: ( 0, 0 ), ( 0, 1,000 ), ( 1,000, 0 ), and ( 1,000, 1,000 ).

import arcpy
import fileinput
import string
import os

from arcpy import env
env.workspace = "D:/UW/Winter_2016/501/w5_more_arcpy8_9_10/ex8"
env.overwriteOutput = True

infile = "D:/UW/Winter_2016/501/w5_more_arcpy8_9_10/ex8/coordinatesEx8.txt"

for line in fileinput.FileInput(infile):
    ID, X, Y = string.split(line, " ")
    array.add(arcpy.Point(X, Y))
cursor.insertRow([arcpy.Polygon(array)])
fileinput.close()

outpath = "D:/UW/Winter_2016/501/w5_more_arcpy8_9_10/ex8"
newfc = "newFeatureClassEx8.shp"
arcpy.CreateFeatureclass_management(outpath, newfc, "Polygon")

cursor = arcpy.da.InsertCursor(newfc, ["SHAPE@"])
array = arcpy.Array()

del cursor

coordinatesEx8.txt looks like this below

1 0 1000
2 1000 0
3 1000 1000
4 0 0

Hanging up on this line:

for line in fileinput.FileInput(infile):

Error message:

ValueError: too many values to unpack

any suggestions? I'm at a loss, I thought it was the set up of my .txt file but keeps hanging on that error no matter what. Not sure how to proceed, this is the question before my actual homework question which will be much harder...also originally the for loop was:

for line in fileinput.input(infile):

which gave a different error...


texts: Coordinates for specific points calculated in QGIS Field Calculator?

texts:
I have got a layer in QGIS that is comprised of connected points which define regions of a country, thus the attribute table consists of 38 regions while the number of nodes in this layer is about 14000. The CRS system is EPSG:3410.

I want to do the following. For each regions calculate the rightmost point and then I want QGIS to give me the coordinates (XY) of this point. I have looked through many forums but I haven't found the answer to this simple problem.

I can find the x-coordinate of the rightmost point, in the field calculator for the regions layer I type x=xmax($geometry), but I do not understand how to get y coordinate of this point.


texts: Proof low deviations in DEM compared to point dataset

texts:
I am not sure how to achieve the following consistently:

I have a regular distributed point dataset (surface elevations). From this I derived several DEM using QGIS IDW Interpolation method. Now I get different results based on the value of the distance coefficient I enter. I examined values between 1.0 and 3.5, where lower values smooth the surface more and larger values increase surface roughness and so terrain features show up more clearly.

In the next step I compared single values from my DEM rasters with the value at the same position in the original point dataset (which I assume to be the "correct" real world values, I want to be represented in the DEM too). I found larger deviations in the rasters with a lower distance coefficient, due to the surface smoothing.

Now I want to ensure, which distance coefficient fits best for my calculation, thus trying to calculate mean deviation between one raster and the point dataset over the entire area.
I though I could just calculate raster - point dataset, then calculating mean value on the resulting dataset, but this wont work with the raster calculator.

I am unfortunately also not too much of an expert in QGIS.

Is there some function achieving the above or a workaround?

Googling, I could not find anything that would do that for me.


texts:
I am using GeoWebCache with GeoServer.
suppose I have seeded a layer to store its tiles. After that around 60% tiles have been created. Now when I requests for the layer from map client, does GeoServer renders the missing tiles by generating these dynamically at it does if tile caching is not used? Or the tiles that have not been created are missed from the layer?


texts: SAD69 conversion into latitude and longitude

texts: I have an electricity distribution database in which the coordinates used are in the SAD69 system. I'd like to convert this data to longitude and latitude, but I do not find formulas anywhere. Could someone explain how I do this conversion, or indicate where I find this information? For example, the coordinate in SAD69 (x = 483473, y = 7415400) would correspond to which coordinate of longitude and latitude?


texts: I am working on a map rendering application and would like to hook it up to a current GIS database to get shapefiles and other data. One possible data resource I have looked at is OpenStreetMap. For my application I would want to have the data locally for processing, so I would need to mirror the database. I know that I can get a dump of the whole database, but ideally I would like to get dynamic updates to my dump as the data is improved; hence I am looking for more like a mirroring functionality. Does OpenStreetMap allow for their whole database to be locally mirrored by users, or does it only allow localized queries?


texts: SQL Spatial select record where intersects another record

texts:
I have two tables - 'Tree' (point) and 'ServiceArea' (polygon) in the GIS database.
Trying to select all Tree records that intersect a ServiceArea where the value in Tree.ServiceAreaNumber field doesn't equal the value in the intersecting ServiceArea.ServiceAreaNumber field.

Something like this:

SELECT
    T.TreeID,
    T.ServiceAreaNumber AS RecordedServiceArea,
    S.ServiceAreaNumber AS ActualServiceArea
FROM
    TREE T,
    SERVICEAREA S
WHERE
    T.STIntersect(S) AND T.ServiceAreaNumber != S.ServiceAreaNumber

The above code is how I see it working when I think about it, I just haven't been able to find any examples of how to actually do a spatial query between two tables.  All the examples I've found have declared a textual point:

DECLARE @g geometry;
DECLARE @h geometry;
SET @g = geometry::STGeomFromText('LINESTRING(0 2, 2 0, 4 2)', 0);
SET @h = geometry::STGeomFromText('POINT(1 1)', 0);
SELECT @g.STIntersects(@h);

EDIT:
Perhaps I need to look something along the lines of

declare @polygon geometry = 'POLYGON((-9486683.581 4810152.256, -9282073.762 4821688.121, -9262037.786 4625578.413, -9477576.319 4628614.167, -9486683.581 4810152.256))'
select * from tablename
where tablename.geom.STIntersects(@polygon) = 1

from SQL Server Spatial Querying however I'm unsure how to specify the second table for the intersect rather than declaring a single feature.


texts: Seeking Open Source ISO Metadata Editor?

texts:
I am looking for a free open source ISO metadata editor for the development of an enterprise database.

It will mostly be used for documenting shapefiles and GeoTIFFs.

Hoping that it will be an intuitive software since a range of skill levels will be working with it.

Any suggestions?


texts: Using spatial index in spatialite query?

texts:
I posted this question originally Using nearest neighbor query in SpatiaLite?

Using the answer given I was able to get the query working on some test data. However I am having trouble applying a spatial index (which is a must as I have over 1 million points and I am trying to find the closest line in a layer that has over 300,000 lines). I've tried the query below on a test data set, but it runs slower than the query statement without the spatial index. What am I doing wrong?

UPDATE UgandaBikeTaxis_test SET dist=(SELECT
    ST_Distance(p.geometry, l.geometry) AND
p.ROWID IN
(SELECT ROWID FROM SpatialIndex WHERE f_table_name= "ugandabiketaxis_test"  AND search_frame=l.Geometry)
    FROM osm_filtered_prj_test AS l, ugandabiketaxis_test AS p
    WHERE p.uniqueID = ugandabiketaxis_test.uniqueID
    ORDER BY ST_Distance(p.geometry, l.geometry) LIMIT 1);


texts: How can I check the resolution of my .img raster file in ArcMap?

texts: I'm just learning ESRI's ArcMap and I just wanted to know where I can find the resolution of a .IMG raster layer.


texts: Get part of polyline based on 2 points in Arc Engine

texts:
Currently I have a layer with polylines; in another layer, I allowed the user to create 2 points on the polyline.

I want to know if there is a function for the polyline feature that I can use to get  the clipped portion of the polyline based on the 2 user-added points. Thanks.


texts: PostGis - ST_AsRaster returning raster smaller than reference raster

texts:
I am trying to convert a point geometry to a raster using a reference raster:

SELECT ST_AsRaster(my_points, ref_raster)
FROM my_table

This works fine except it is only returning a subset of the reference raster:

SELECT (ST_Count(my_rast))::numeric as pixel_count,
       (ST_Count(ref_rast))::numeric as ref_pixel_count
FROM (SELECT ST_AsRaster(my_points, ref_raster) as my_rast, ref_raster
      FROM my_table) AS my_inner_query

returns a pixel count of 2739 for my created raster and 3444 for my reference raster.

I have tried to understand the documentation but it is not very clear to me. Does anyone know the correct way to create a raster from a geometry using a reference raster so that the created raster is the same dimensions?


texts:
I have two sets of MODIS imagery, changed and re-projected to GEOTIFF:
1.includes the south part of my region (100 images)
2.includes the north part of my region (100 images)
I want to merge them into one file for each date.
for example at May 1st, I have one image covering north, one covering south...and I want to merge them into one... Also for May 2nd, 3rd and for all later dates.
Any idea how to do that automatically?


texts: Clipping virtual raster generates out of memory error

texts:
I'm using QGIS 2.16.1 on Windows to try to extract elevation data into curves that I can use to cut/build a wooden topo map.  I'm vaguely following this guidance, but also am changing things up a bit as needed (my data is across multiple files, etc).  I'm in really deep water here as I'm a woodworker and not a mapping person (yet!).

I've imported 6 raster layers of elevation data from the ArcGRID datasets, merged them into a Virtual Raster, and ran Raster > Extraction > Contour at an interval of 125.  After several minutes of chug, I get the contours.  Awesome!

My problem: Trying to save myself some data chunking time (I'm doing some smoothing and eliminating of too-small-for-wood contours), I wanted to eliminate the data and contours that are outside the boundary of the region I wanted to topo (several of the data files only overlap my region by a tiny amount).  I loaded a shapefile for the region boundary (it's a national park so finding that was easy) and used Raster > Extraction > Clipper to clip the Virtual Raster to the shape.  It appears to work and generate a new Clipped Raster (my name for it) in the right shape that visually matches the original data.  Problem is, when I try to run Contour on the Clipped Raster, it generates an ERROR 1: CPLRealloc(): Out of memory allocating 416 bytes.

As first I assumed the clipping shape was too complex, so I tried with just a square and it fails with the same error.  Seems any Clipping of the Virtual Raster drives the error.

Is there any way I can fix this or any other method to take the elevation data from multiple files, "combine it" and extract contours?


texts: Mapbox Natural earth map disappears on zooms 0 & 1

texts:
As the title says, I'm currently using mapbox to create a little atlas app the issue is whenever I zoom out to anything lower than zoom level 2, the map dissappears with the map background being the only thing visible. Here's the CCSS I'm using

Map {
  background-color: #F2F7F7  ;
}

#ne_10m_admin_0_countries  {
  polygon-fill: #E3E6E8 ;
  line-width: 1;
  line-color: #C2C6C6;
}

Any ideas why this is happening?


texts:
I'd like to read data (from a device connected to the serial port) in QGIS using Python. Its not a GPS device; it is an AIS receiver so the QgsGPS tools obviously won't work.

I have got as far as;

QgsGPSDetector.availablePorts()

Which does appear to recognise the device;

[(u'localhost:2947:', u'local gpsd'), (u'\\\\.\\COM11', u'COM11:')]

I understand that serial/pyserial library is not included in the PyQGIS environment. If it was I think I could figure it out easily.

Any suggestions how to write a short script to read raw data from the serial/COM ports?


texts: How to derive vehicle paths from timestamped location data?

texts:
I've got a bunch of location data for snow plows. There's a unique id for every plow, a timestamp, and x & y coordinates.

Instead of knowing where hundreds street plows were at one minute intervals, I'd like to know which streets have been plowed.

I know there are ways of turning points to polylines, but does anybody know of a method to snap the points to the street grid? In other words, of deriving vehicle paths that hug the streets?

The reason I think that turning points to polylines probably isn't enough is because I only have location data for every plow every 60 seconds. This means there aren't coordinate points for when plows were turning at certain intersection, so if you just drew lines between points some of them would cut through blocks.


texts: Using point shapefile as training sample in GEE

texts:
I want to do random forest classification. I prepared a point shape file with training sample outside GEE, the true land cover class being in the field code. Then I imported the shapefile in GEE (as training_shape) and wanted to sample a raster image (input_img_m) in order to get the raster values for each point and thus generate training sample for classification.

However, I am not sure which tool to use. I tried with .sampleRegions. I know this is working exactly I want it with polygons and I hoped it will work with point as well. But after 20 min of processing I get error: Error: Unable to export features with null geometry.

This is part of my script:

var training = input_img_m.sampleRegions({
  // Get the sample from the point FeatureCollection.
  collection: training_shape,
  // We'll classify on 'code"
  properties: ['code'],
  // Set the scale to get Sentinel-2 pixels.
  scale: 10,
  tileScale: 4
});

Export.table.toAsset({
  collection: training,
  description:'exportToTableAsset',
  assetId: 'samples',
});

My question is what is the better way to use training sample shape file for classification in GEE if the geometry of the shapes is point?


texts: Can I change the colors in FME translation log?

texts:
Due to my small colorblindness I face difficulties in recognizing red text in black text.

In the FME Translation log errors are always shown in red. Is it possible to change this color (or that of the black text)? Because in this setting I can't use the software. I couldn't find it in the FME Appearance options, perhaps someone has a good solution.


texts: Python code help for ArcGIS Field Calculator tool (Model Builder)

texts:
Is there any python code that:

1st, sort a column (example name by alphabet)

2nd, grab the first cell of sorted column and assigned it to a variable like A=..

I need all of this to be done in field calculator and inside show block. Please give me a sample or help.


texts:
The FME transformer "NeighborFinder" seems to have no problems finding the nearest point to a polygon. But I want it the other way around. I want to set the points as the "base" and the polygons as "candidates". But as a result I only get the same coordinates of the points.

How can I find the polygons nearest to the points?


texts: Identifying top 10% of cells from raster for number of regions using ArcMap

texts:
I have a raster with different regions (regions in the geographical sense) and I would like to create an output, that shows me the top 10% of cells for each region. I found an answer on how to do it for a single raster, but not for multiple regions within one raster (Identifying top 10% of raster cells using ArcGIS Desktop).
Has anyone got an idea on how to do it or do I need to split up the raster first and then run the calculation for each subraster?
I am using ArcMap (Spatial Analyst included). The raster is a representation of agricultural soil suitability of a state. I have a seperate feature class with 25 subregions of that state, for which i would like to identify the top 10% of best suited cells for each subregion (10% of the agricultural area of each region with the highest suitability). The image shows the raster (red to blue is the suitability, white is no data), black lines are the subregions of the state.


texts: How to write an equivalent to MapBasic's CreateLine in pyQGIS?

texts:
I know this kind of question has been asked, in part, before - specifically: Seeing QGIS Python Commands and: How is object created as line using createline stored in .tab file?
but this is a specific plea for help!

I've lost my (ancient) version of MapInfo, and so am attempting to recreate some analysis in QGIS. I relied heavily on the following MapBasic Command:

set coordsys table "yourtable"
update "yourtable" set obj = createline(x1, y1, x2, y2)

It joined every single point to every single point in the table and retained attribute data, thus giving me the ability to subsequently analyse the file. I loved this command for its' sheer simplicity.

I'm currently attempting to teach myself some Python, but from a very low knowledge base - so, I was wondering if anyone had seen something similar scripted in Python?!?

UPDATE: Thanks to Nathan and gene, comments really helpful!


texts: How to display SQL view query (from GeoServer) in OpenLayers

texts:
I have written this SQL view code in GeoServer:

select *
from pilot
where ((Price between %MinPrice% and %MaxPrice%)
and (area between %MinArea% and %MaxArea%));

So according to the changed parameters this Request is produced as follows:

http://localhost:8080/geoserver/pilot/wms?service=WMS&version=1.1.0&request=GetMap&layers=pilot:FilterAreaPrice&styles=&bbox=535410.289995867,3949726.85999973,536327.370000143,3950801.26000122&width=655&height=768&srs=EPSG:32639&viewparams=MinPrice:1;MaxPrice:3;MinArea:1000;MaxArea:3000&format=application/openlayers&

My goal is that the user enters parameters in a form and then clicks on "submit" button and then results are displayed based on the parameters entered in my OpenLayers interface.

.

Additional information is attached


texts:
I´d like to construct a shapefile from a Pandas Data Frame using the lon & lat rows.

I have got a csv file and I process it with pandas to make a data frame which is easier to handle

Is it posible to do that without make a loop line by line ?


texts: Filtering 2 tables that are in relations (1:n)

texts:
I work with two tables which are in relations from 1 to n (one to many) : An individual owns several goods.
My individuals are in the main form and their assets are displayed as a table in the sub form.
I want to map the doctors who own a pet for example.
How do I do the filtering?
For now, I make a query on each table and then I join them by an attribute join. It works well.
Is this the right way ?
Can we do this more straightforward?


texts:
I'm developing a map in an Android WebView, using OpenLayers and a WMS Service. I can only zoom a certain amount of times (which I would like to increase).

My code is:

map = new OpenLayers.Map({div:"map"});
baseLayer = new OpenLayers.Layer.WMS("WMS", link, {layers:"basic"});
map.addLayer(baseLayer);

I've tried from the page Setting Zoom Levels

var options = { numZoomLevels: 25 };
map = new OpenLayers.Map( {div:"map"} , options);

but then the map doesn't even show. And this:

map = new OpenLayers.Map({div:"map"});
baseLayer = new OpenLayers.Layer.WMS("WMS", link, {layers:"basic", numZoomLevels: 25});
map.addLayer(baseLayer);

but in this case, nothing changes (I can only zoom in 16 times, as before). Where is my mistake?


texts: Postgis and python parameterised queries

texts:
I have a python script that i am trying to run in order to dump a postgres table and I noticed there is a pgsql2shp command in bash to perform this easily but I would like to incoperate the same command in a python script and also use parameter substitution. I thought the following would work but it hasn't yielded any result.  What is the equivalent of the this command in python

(pgsql2shp -f %s.shp -h localhost -p 5432 testdb "%s_%s")%(table,table,test)


texts: How much Water on Earth(USGS) : A Wrong Map representation?

texts:
I happened to read a USGS article, How much water is there on, in, and above the Earth? showing a representation of all the water(which includes water from Oceans, Seas,Bays,Ice caps, Glaciers, Permanent Snow and ground water) displayed as a small sphere as shown below.

I find it hard to believe although the site has given good statistics for it. Is this a correct representation or is it a map gone wrong?


texts: Getting from a USGS historical topographic map to a Unity terrain heightmap

texts:
Is it possible to start from a USGS historical topographic map and get to a Unity terrain heightmap?  To see an example of the kind of historical topo map I am talking about:

Go to: https://ngmdb.usgs.gov/topoview/viewer/#9/35.8690/-109.0269
Enter "Fort Defiance, AZ" in the upper right search box (with Location tab selected)
Slide the right date slider down to "1886"
Zoom out in the main map to see the selected area

If you then click on the image of the map in the right panel, a preview will come up showing the 1886 topographic map.  The map can be downloaded as a GeoTIFF, which I did.

I have next to no experience with GIS, but I've tried a number of things with gdal, QGIS, Mapbox Studio, among others, based upon steps others have taken, but to no avail. I started with attempting to follow this:

Importing DEM Terrain Heightmaps for Unity using GDAL

then this:

Cropping GeoTIFF to square at specific lat/long for use as Unity 3d terrain heightmap?

I cannot get anything to work. I'm thinking that, while the GeoTIFFs of these historical topo maps may contain some lat/long information, there is no elevation data and, in that respect, it is simply an image.  Am I correct in this?

Here is the gdalinfo output for the example map above:

~/workspace/unity/terrainmaps/ -> gdalinfo az_fd_geo.tif
Driver: GTiff/GeoTIFF
Files: az_fd_geo.tif
Size is 4958, 5954
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["NAD27",
        DATUM["North_American_Datum_1927",
            SPHEROID["Clarke 1866",6378206.4,294.9786982138982,
                AUTHORITY["EPSG","7008"]],
            AUTHORITY["EPSG","6267"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4267"]],
    PROJECTION["Polyconic"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-109.5],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-52090.740425979456631,3992327.326727173291147)
Pixel Size = (21.166666666666664,-21.166666666666615)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=300
  TIFFTAG_YRESOLUTION=300
Image Structure Metadata:
  COMPRESSION=YCbCr JPEG
  INTERLEAVE=PIXEL
  SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left  (  -52090.740, 3992327.327) (110d 4'41.43"W, 36d 3'41.62"N)
Lower Left  (  -52090.740, 3866300.993) (110d 4'12.34"W, 34d55'32.52"N)
Upper Right (   52853.593, 3992327.327) (108d54'48.09"W, 36d 3'41.47"N)
Lower Right (   52853.593, 3866300.993) (108d55'17.61"W, 34d55'32.38"N)
Center      (     381.426, 3929314.160) (109d29'44.87"W, 35d29'42.08"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  Overviews: 2479x2977, 1240x1489, 620x745, 310x373, 155x187, 78x94, 39x47
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  Overviews: 2479x2977, 1240x1489, 620x745, 310x373, 155x187, 78x94, 39x47
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  Overviews: 2479x2977, 1240x1489, 620x745, 310x373, 155x187, 78x94, 39x47
~/workspace/unity/terrainmaps/ ->

So, is there a way to do what I'm trying to do, perhaps a different approached based on whatever the file has to offer?


texts: How to divide line into three equal parts?

texts:
I have a problem with line division. I'm writing a program in JS/C++ generating KML files.

I was trying to divide a line using my own function (which is correct on cartesian space):

    var part1Lon = startLon+((-startLon + endLon)/3);
    var part1Lat = startLat +((-startLat + endLat)/3);

But when i visualized point [part1Lon, part1Lat] on Google Earth - and it's not on the line! (As I could have expected)

So, second attempt was to convert longitude and latitude to Mercator first, do calculations, and convert coordinates back to WGS. It produces similar result - point is not on the line.

What is correct algorithm or the correct projection for doing such operations?


texts: Calculating minimum distance between all points of 2 different point layers using QGIS?

texts:
I use QGIS 2.8.3 or 3.4 whatever works better.

I need to calculate the mimumum distance between all points of 2 different point layers. I was reading this post https://gis.stackexchange.com/a/62212/147194 and used the distance matrix tool accordingly. Now I am unsure what the output layer contains. The point input layer had 488 feautres and the target point layer had 605 features.

But why does the output layer contain of 3801 ?

The last column of the output layer probably gives the distance but in which measure unit?

I need to 2 this process for 2 other point layers and want to be sure that I am on the right track before I continue.

Is 'Distance Matrix' the easiest way to calculate the minimum distance between all points of 2 different point layers?

And what does the number at 'use only nearest (k) target point' stand for?

Following the instructions (see link above) I put '1' but without knowing why?


texts:
it's the first time I'm working with MODIS data and I am a bit confused about how bits work and how masks are created. After a bit of research I've been able to understand, more or less, how bit filtering works -- using the help found here and MODIS documentation files (1 and 2).

I am using MOD09GA product to create cloud and aerosol free mask data. I built my binary filter based on the information found in the next table (file 2, page 24).

When doing bitwise operations in python and applying the bit filter to the QA band (state_1km), I am not fully sure that what it returns is right, since it generates a mask with 4 integer values (0, 1, 2 and 3) when I was expecting 2 (0 and 1).

I will show this with an example. You can download the .hdf file from here, although you'll need you login credentials.

import gdal
import numpy as np

# Opening MODIS hdf file
 modisImage = gdal.Open("data/MOD09GAMOD09GA.A2008001.h10v05.006.2015169041625.hdf")

# Getting subdatasets
subData = modisImage.GetSubDatasets()
data = {}
for fname, name in subData:
    newname = name.split("] ")[1]
    without_type = newname.split(" (")[0]
    print(without_type)
    data[without_type] = gdal.Open(fname).ReadAsArray()

# =============================================================================
# # num_observations_1km MODIS_Grid_1km_2D
# # state_1km_1 MODIS_Grid_1km_2D    <-- I want this one
# # SensorZenith_1 MODIS_Grid_1km_2D
# # .
# # .
# # .
# # .
# # q_scan_1 MODIS_Grid_500m_2D
# =============================================================================

# now call the desired band
qc = data["state_1km_1 MODIS_Grid_1km_2D"]

Cloud state, cloud shadow and land/water bit masks are returning a mask with the values expected (0 and 1), but when I apply the same approach to the aerosol quantity I can't understand the mask returned.

# Aerosol quantity information is in the Bitshift 6 and Bitmask 0-3 (first and third columns of the table). In this case, I want to return only high quantity aerosol

qcAerosol = (qc >> 6) & 3

# Checking values returned (more than 2 values)
np.unique(qcAerosol) # array([  0,  1, 2, 3], dtype=uint16)

# Taking a closer look to te values returned
bin(1)  # '0b1'
bin(2) # '0b10'
bin(3) # '0b11'

The qcAerosol mask created is the one in the image. As you can see, it has four colors corresponding to each value.

Each value returned in the mask, maps a binary value in the table:

1 corresponds to the bit mask for 'low aerosol quantity'
2 corresponds to the bit mask for 'average aerosol quantity'
3 corresponds to the bit mask for 'high aerosol quantity'

Since I want a 'high aerosol quantity' mask, I need a last step to create one.

qcAerosol[qcAerosol != 3] = 0
qcAerosol[qcAerosol == 3] = 1

Getting the next mask.

Python qcAerosol

However, as I noted above, bitmasks for cloud state, cloud shadow and land/water only return two values, 0 and 1. I think there could be a problem with my approach, therefore my question is raised: When working with MOD09GA product, is supposed QC bitmask for high aerosol quantity to return 4 values?


texts: How to embed a coordinate pair within URL to search OpenStreetMap, Google, Yahoo or Bing Maps

texts:
Is there a way to embed a latitude longitude coordinate pair within a url link to search for a location on either OpenStreetMap, Google and/or Bing Maps?  I'd like to use the url link within Excel, LibreOffice and within an html pages.  I have tables with coordinates that I'd like to configure with a link to show its location on a map, without building an online mapping application.

If not, then is there a simple javascript that I can encapsulate the coordinate within and fetch that from OpenStreetMap, Google and/or Bing Maps?

It seems like the url would be something like the following, but I can't get this to work:
http://www.openstreetmap.org/?30.486,-90.195


texts: Why is ArcGIS Pro so much faster than ArcMap, especially when dealing with raster?

texts:
I began using ArcGIS Pro back in April or so. First I didn't really like the ribbon layout as it gives me a strong vibe of Microsoft products. But learning new things has always been exciting for me, therefore I kept at it.

Over the last couple of months, I found that ArcGIS Pro is so much faster at handling raster. For example, when I used ArcMap to do Extract by Mask on high-resolution aerial files, I had to leave the computer on after work so that it could be finished the next day. Sometimes I even had to leave it on for the entire weekend. However, now with ArcGIS Pro, similar tasks won't take longer than 10 mins. That difference is insane to me.

Does anyone else have similar experience? I'm wondering if anyone knows what fundamental changes Esri has made in ArcGIS Pro. Did they change the algorithms for most tasks? I find it unlikely since the ArcPy code is the same.

I don't have any software development experience, but I really want to learn a little bit more about it. Can anyone share with me how such a dramatic change transpired?


texts: I am working on a plugin in qgis and i have to make a document from A to Z about that but on French . and I am looking now for some documentation of python,QT , pyqt4, and QGIS in French.


texts: Understanding if pixel value is brightness or height?

texts:
I have geotiff. How to understand if pixel value is brightness or height?

I am downloading SRTM data from
https://earthexplorer.usgs.gov/

And trying to open them on ArcScene, but when I am trying to make them look as 3d I am getting some crap, like pixels do not have height, but just color value.

My image become look like column.


texts: Dojo Slider to control tiled map services?

texts:
am creating a web map to show Sea Level Rise in my County. I want to use a Slider for a user select the increment of Sea Level Change that wish to see.

I have been able to create Programatically create the Slider using:

var vertSlider = new VerticalSlider({
       minimum: 0,
       maximum: 10,
       discreteValues: 1,
       //value: 0,
       intermediateChanges: true,
       style: "height: 200px; width: 15px;",

       onChange: function (value) {
           dom.byId("sliderValue").value = value;
       }
       }, "vertSlider");

The maps layers are created with:

var layer1 = new ArcGISTiledMapServiceLayer("http://url");
var layer2 = new ArcGISTiledMapServiceLayer("http://url");
... and so on till 10 layers are added.

I need now to be able to pass the value of the slider to draw one of my map layers.  So when the slider is moved from position 1 to 2, layer 1 would be turned off and 2 would display.

How would i do this using JS or the ESRI API?   There are many examples using the TimeSlider, but those do not translate to the Dojo Slider.


texts: Saved maps in geoexplorer (geoserver layer security)

texts:
I am a bit confused with the save maps function in geoexplorer.
I have a vps server with the opengeo suite 3.1, I added a table to the vps local postgis db server to the geoserver.public database. I registered this table as layer in geoserver and it's working fine in geoexplorer composer. My problem is when I save the map in geoexplorer and log out and log in back (as administrator) I cannot use the query button and the edit button (they are inactive). It is working fine (so I can use the query button on saved maps) if I don't change the layer security in geoserver. But I changed because I don't want to allow to anonymous users to see any layers.

with this settings: the anonymous user cannot see the layers preview in geoserver and in geoexplorer. Admin can see every layer after logged in. But If I save the map (not publish), the query button and the edit button will be inactive.
..r   ADMIN,GROUP_ADMIN
..w   ADMIN,GROUP_ADMIN

with the default settings:
..r   *
..w   *
The anonymous user can see the layers, admin can use the query button on saved maps.

I would like to hide the layers from the anonymous users,
and I would like to use the query button on saved maps.
Is it possible to do both of them at the same time? I'm sure is possible, just I don't know how can I do that.

I cannot solve this problem and it's very frustrating. It would be great if somebody can give me some idea to solve this issue.

Thanks,
encleadus


texts: Merge 2 layers based on an attribute column in QGIS

texts:
I have a vector shapefile map that is divided into polygons. I want to add two columns to the attribute table with the Eastings/Northings coordinates of each polygon centroid.

So far what I've tried is:

menu > Vector > geom tools > polygon centroids

This gives me the centroids, but as a new layer. What I can't figure out how to do is to merge this layer with the original.

There is a command: menu > vector > join attribute by location.

But what I want is, I guess, 'join by attribute column' or something.


texts: How do I filter unique values using a definition query on a File Geodatabase?

texts:
Is there a way to filter unique values only from a column in ArcMap using definition query?  i.e to use a Definition Query to filter out just those values which occur only once.

In a Personal Geodatabase I can use "DISTINCT", what about File Geodatabase format?

Thanks for your help.


texts: Changing precision on a column in the attribute table without losing data?

texts:
I have data in an Excel file that I want to join with a shapefile in my project.
The Excel file contains values like volumes. The Excel file and the shapefile I want to join it with both contain the same property names.

The problem is, when I join the two, I get too many decimals in the attribute table. I get six, but I only want it to be one.

I tested to change the format of the cells to one decimal, and then I saved the file as *csv. But, somehow I seem to loose data in the process? It should be 42 rows that contains data in the column i am working with, but I end up with only like 10?

Another thing is that my property names contains the letters Å. Ä and Ö, perhaps this complicates things?

I read the question How can I change the precision in the attribute table?, but it doesn't mention the problem with the data loss.

My question was tagged with qgis, but I didn't write what version I'm using. I use QGIS 2.2.0. Would it help to upgrade to the latest version?


texts: Remove unused nodes from OSM data files?

texts:
For some experiments on OSM data size, I am using osmfilter to remove ways and relations with certain tags from the planet dump. But AFAIK osmfilter does not remove the corresponding nodes. So I'm left with a huge number of nodes that are not used by anything, and that distort my data size results.

Is there a tool to remove unused nodes (those that have no tags and are not referenced by either a way or a relation) from OSM data files (e.g. .pbf or .osm.bz2 files)?


texts: I have created 2 Categories in a feature using classification.Now how do I convert them into new shapefiles?


texts: How to change font used by ArcGIS Pro controls and dialogues to Khmer?

texts:
In Cambodia, we need to use Khmer font to show our own language characters. ArcGIS Pro which has Ribbon interface support the display of the Khmer font but we need to change the font to use Khmer font. I am wondering if it is possible to change font used by Ribbon controls and dialogues to the one that support our Khmer Unicode?

In the picture, Khmer displayed properly when we change font to Khmer Unicode while the application dialogues show rectangles, sign of unsupported font.


texts: Openlayers - Problem with selecting features from overlapping vector layers

texts:
I face the following problem.
I have two vector layers. One has linestring and the other point geometry.

When I add this line to my point geometry vector layer:

          pointLayer.setZIndex( 1001 );

I get the point layer on top and then I am unable to select any feature from the linestring vector layer.

How can I solve this problem so all the features are clickable? I wouldn't expect this behavior as the point vector layer just has some points.


texts: Metadata edit bar disappeared in Arc Catalog

texts:
could anybody help me? in my Arc Catalog the edit bar for Metadata disappeared.
In my Catalog it looks like this:

But this ist from Esri website, it should look like this:

How can I add this Editing bar?


texts: Combine a checkbox with a Split panel in Google Earth Engine to compare two images

texts:
I was wondering is someone can tell me if it's possible to include a checkbox within a Split Panel. In my case, I'm comparing data from Sentinel 5 with nightlights. The default option is the images of NO2 from the Sentinel 5. The right panel shows data from 2020 and the left panel from 2019. When I select the checkbox, I would like to see the images from nightlights for the same periods, instead of NO2. If unchecked, it shows NO2 data.

Here is my code: https://code.earthengine.google.com/a84bf389db81358f5c42a5a422ddfe71
I'm new with this App. Any help will be highly appreciated.
Thanks!