texts:
I have a map with many different layers all lining up fine and projected the same, but a new layer I've recently added is not lining up exactly. It's not off by a lot but you can definitely see the boundary and the roads not aligned. If they are in the same projection, shouldn't they be overlapping perfectly? Or is it ok to have them slightly off? Pics for reference below:
Data source for both layers in the picture:
Projected Coordinate System:    NAD_1983_UTM_Zone_18N
Projection: Transverse_Mercator
False_Easting:  500000.00000000
False_Northing: 0.00000000
Central_Meridian:   -75.00000000
Scale_Factor:   0.99960000
Latitude_Of_Origin: 0.00000000
Linear Unit:    Meter

Geographic Coordinate System:   GCS_North_American_1983
Datum:  D_North_American_1983
Prime Meridian:     Greenwich
Angular Unit:   Degree


texts:
This is a bit of a beginner question I guess but I just can't find a clear answer anywhere...

I have a projected satellite image (warped via gdal) and would like to set the no data pixels to a transparent value so that when the wms layer of my (image) coverage is added on my map only the projected image is displayed (and not the balck area around).

Is it in the sld for my layer? is it in the wms parameters? do I need a rgba image format? etc... Any pointers would be highly appreciated!


texts: Programmatically download NED, Landcover, etc. data from USGS

texts: What are my current options for downloading geodata (NED, landcover, etc.) from the USGS programmatically?


texts: How can one get the value of the point where the AxMapControl is centered at?

texts: I am using the AxMapControl to display a map and I want to record a point where I am and want to be able to return to the same point later on using the CenterAt() method.


texts: Making and exporting high resolution (detailed) map in QGIS

texts:
My co-worker gives me a map in JPEG format which is exported from ArcGIS Desktop. When I zoomed it, it shows every detail of the map.
I am using QGIS and I export a map from QGIS. When I try to zoom it, I can't see the detail of the map like the ArcGIS one.
How can I make it visible in detail way like the ArcGIS map?

QGIS Resolution with 2 point font size.

ArcGIS Resolution.


texts: I am facing problem merging vector layers (point layer,and buffer layer)


texts: Grouping large dataset of lakes into n clusters based on target total perimeter?

texts:
There's a longer explanation of why I'm interested in this below, but in short, I'm wondering how to run something like a k nearest neighbors formula, with the number of groups determined by a maximum total characteristic(perimeter) in each group.

This is part of a cutesy project a friend and I are working on-

We're from Minnesota, often called the land of 10,000 lakes and are exploring what it'd look like to actually run around 10,000 lakes over a long period of time.

Finding the 10,000 largest was pretty trivial, as was generating perimeter measurements.

To start making "courses" to run, we're looking at dividing the lakes into groupings that are:

1) close together
and
2) about 8-15 miles in perimeter total between the lakes.

Any ideas on how to best accomplish this? I'm playing around with grouping analysis in ArcMap right now, but I'm pretty sure that there's a more optimal solution someone's made already.


texts: Can we use PostGIS functions in Openlayers?

texts:
i have a traffic data (vehicles gps location(Point)) for every minute in a specific location

sample data : Vehicle_No , Long, Lat

and i need to show the real time traffic for that specific location in openlayers  which has to be updated every minute.

one solution is to import traffic data into PostGis and create a new cloumn which use ST_Buffer function to find out how many no. of vehicles are near by each other and can export the table to the openlayers and give them colour (red,orange,blue) according to no. of vehicles (the buffer column produced).
- the problem with this approach is that we have to manually do the buffer operations and update the openlayers every minute.

Can someone suggest me some better approach to do this?


texts: Is there cartographic convention to ordering legend items (high-low or low-high)?

texts:
When you symbolise a layer in ArcMap, the symbols are ordered from low value to high value by default. Is that necessarily "correct"?

I've never officially studied cartography - is there a convention on whether to order items from low-to-high or high-to-low in a legend? Are both ways correct?

Edit: here's a screenshot of the symbology dialog in ArcMap, showing median weekly rents. It's ordered from low-to-high by default, but has the option to flip the order:

Here's how it appears via ArcGIS Server on a website:

Note that my question is about cartographic standards in general, and is not limited to ArcMap.


texts: Where to put vector layers with different CRS in a PostGIS database?

texts:
We have a PostGIS database with vector data with CRS British National Grid EPSG:27700. It is accessed remotely by users with QGIS.

I now need to set up a similar vector layer for Northern Ireland, which has CRS Irish Transverse Mercator EPSG: 2157

What things do I need to consider when I am deciding where to put a vector layer with a different CRS in a PostGIS database?

I have had a look at this post relating to raster data, but its very detailed and specific to that particular question How should I manage PostGIS Raster data with different projections?


texts: Creating square of exact size (e.g. 1 km²) around point in Google Earth Engine

texts:
I struggle to create a geometry with the exact size of a square with for example 1km² around a specific GPS point, to later export an image.

Is there a way to define the length of each side of the square and center it around a given location?

Or if not, can I calculate the corners of the square by using my GPS data and create the square like this?

// Create a geodesic polygon
var geometry = ee.Geometry.Polygon
([ [[-105.35, 39.95], [-105.35, 40.05],[-105.2, 40.05], [-105.2, 39.95], [-105.2, 39.95]] ]);

// Display the polygon on the map
Map.centerObject(geometry);
Map.addLayer(geometry, {color: 'FF0000'}, 'geodesic polygon');

// Print polygon area in square kilometers.
print('Polygon area: ', geometry.area().divide(1000 * 1000));

This is what I can do. Basically just creating an Polygon without the same edge lengths (not a square) and measuring the size of it in km². But I am not able to define the exact size of it and to center it over a specific point (sampling site).


texts: Geometry issues when creating table from existing spatial table

texts:
In short, I want to split an existing spatial table into two parts. The split itself is straight forward. I start with this (fully functioning) table:

CREATE TABLE table
(
  my_identifier smallint,
  geom geometry,
  CONSTRAINT enforce_dims_geom CHECK (st_ndims(geom) = 2),
  CONSTRAINT enforce_srid_geom CHECK (st_srid(geom) = 4326)
)
WITH (
  OIDS=TRUE
);

and create the new tables by:

CREATE TABLE public.table_part1 AS
SELECT * FROM public.table WHERE idmy_identifier = 1;

CREATE TABLE public.table_part2 AS
SELECT * FROM public.table WHERE my_identifier = 2;

However, the problem is that although the geom column gets transferred to my new tables, the geometry (?) itself doesn't get transferred correctly:

CREATE TABLE table_part1
(
  my_identifier smallint,
  geom geometry
)
WITH (
  OIDS=FALSE
);

I tried running SELECT Probe_Geometry_Columns(); but it didn't help. Do I need to insert data into the geometry_columns table? If so, exactly what?


texts: Hillshade not displaying correctly in ArcGIS Pro?

texts:
I have a large DEM of England that I have clipped to my study area and have been trying to generate a hillshade from. I tried a few variations with the parameters in altitude and azimuth (leaving z =1), but I continue to get this really grainy result using ArcGIS Pro.

The cell size is 0.00027777 decimal degrees. and I tried to re-sample to 1 but got an error that the size was too large. I tried running calculate statistics and building the pyramids on the original DEM but the result hasn't changed any.


texts: How to find geographic coordinates of intersection point of two lines originating from two different geographic coordinates?

texts:
I have two planes emanating from two different geographic coordinates as per wgs 84 reference ellipsoid namely (lat1, long1, alt1) and (lat2, long2, alt3) with a bearing b1 and b2 with respect to north. These two planes intersect to form a line. This line intersects wgs84 reference ellipsoid at a point latitude, longitude and altitude (lat3, long3, alt3) .
How do I find this intersection point ?
Can some one explain the math behind this ?


texts: Making qgis2web web map display description rather than value of field with valuemap widget

texts: I am making a web map with qgis2web which I want the pop-up to display the description rather than the value from a field which I have a value map widget on. Is this possible, or should I just make another field for my description?


texts: PostGIS: How to find routes near to some point

texts:
So, I have groups of polylines returned from google map API (Green, Red and Yellow), some starting point p (blue one), and a radius r. I need to find all routes whose distance from p is no greater than r (in this case red and yellow)

I know that I can calculate it manually (find equation of each line and compare it to the equation of the circle) but it can be inefficient for many routes. So my question is:

Is there builtin function in PostGIS which will help me accomplish this task?


texts: How to fix unreachable nodes in QGIS GRASS network?

texts:
I have created a network with the GRASS module v.net. But when I try to create the travelling salesman path with v.net.salesman, I get:

  Destination node 9383 is unreachable from node 9398

If I put the topo lines or points from the network into QGIS I can find the two nodes, but I can't manually figure out why one is unreachable from the other. The Topology Checker shows no errors. I know that the v.net.salesman works because I was able to create the path with another dataset. This dataset is faulty, even though v.net created a network. Is there a tool either in GRASS or QGIS that I can use to fix this problem?


texts: Getting attributes from intersecting lines to label generated intersection point?

texts:
I have an urban road network (as a feature class), consisting of only polylines. I have to make a database of all the intersections between the streets.

As I generate those intersection points, I don't know how to address the names of the streets as an attribute for each point. For example, if I have streets with names "1" and "2", the point "x" at their intersection should be named "1,2" or so, if I'm not mistaken.


texts: QGIS 3.10 layout window: Where is the button to filter the legend by map content?

texts: QGIS 3.10: I'm missing the filter-legend-by-map-content-button below the legend in the layout window. How to filter the legend without this button?


texts: How to rollback the committed changes in the Vector Layer in qgis

texts:
I selected few features in the vector layer and updated the field. But I want to rollback the changes committed to the vector layer. How can I do that in QGIS?

I have used the following coed:

layer.rollBack()

But it is of no use. How can I rollback the committed changes in the vector layer in QGIS?


texts: How to perform relational queries against an ArcSDE geodatabase using the ArcGIS API for JavaScript?

texts:
I have a geospatial database which contains both spatial and not-spatial tables
I would like to publish my data on ArcGis for server and perform both simple and spatial queries.
Also, foreign key relationships should be used so as to perform more complicated queries.

To be more specific:
There is a need to be able to execute application specific queries over multiple relational tables and some of them should contain geometries. So, for examle through the arcgis I have to be able to ask a question like : Give all the stations (stations Table) that are connected to ports(ports table) and have color blue (color table). The criteria should be added dynamically so, for example one could select through the jS api, the stations to have an additional characteristic from another table

I know that  you could build the the Query object dynamically, but it's only for attributes of a single table (or joined table). If you could not know in advance which exactly tables you would join (for example implementing a filter) or how many of them you would join (table A with B , A with B and C , A C , A B C E, A E etc).I suppose you dont have to  precompute all the possible join combinations.

Should I publish the entire geodatabase in some way?
Would relates help in that case?
Since queries should be constructed dynamically, I suppose it is not possible to be based on precomputed joined tables.

I need to do this maybe through the JS APi
I am using arcgis for Deskotop 10.2 and oracle DB 11g


texts:
I am trying to include a jpg to a geoext popup.

There are various shapefiles, but only one needs to load pictures. With only openlayers it was easy to do, I just used the content.ftl from geoserver 2.4. However, we switched to geoext and sadly I am not able to figure out how to do it now.

Can somebody give me a hint how to force GeoExt.Popup to use the different content.ftl of geoserver - or is there another way?


texts: OpenLayers - order of layering from a KML file

texts:
I have a set of KML files hosted on a webserver, which two applications use; one uses OpenLayers and one uses Google Maps.

The problem is the order they render the KML. Google Maps (and Google Earth) renders the first linestring in the KML as the top layer. However, OpenLayers places the first one underneath, with later linestrings appearing on top of it.

Hence the same KML file looks different depending on which application is displaying it.

Is there anyway to choose which object appears on top in OpenLayers?

The code which adds the KML to the map is:

function routeMap_addKML (label, kmlurl, show) {

var layer = new OpenLayers.Layer.Vector(
        label,
        {
          strategies: [new OpenLayers.Strategy.Fixed()],
          protocol: new OpenLayers.Protocol.HTTP({
            url: kmlurl,
            format: new OpenLayers.Format.KML({
              extractStyles: true,
              extractAttributes: true,
              maxDepth: 2
            })
          }),
          'onFeatureInsert': function () { reScaleMap(routeMap); }
        }
      )

      layer.setVisibility(show);
      routeMap.addLayer(layer);

return layer;

}


texts: map.on('click') limit number of pins added to map - leaflet js

texts:
I'm using Leaflet js and I want a user to be able to add a single point to the map (and not get trigger happy and add lots of pins), what do I need to do to modify the below so that this can happen.

map.on('click', function(e){
   var Marker = L.marker(e.latlng,{
       draggable: true
   }).addTo(map);

Si


texts: Python new packages installation with ArcGIS interpreter?

texts:
I would like to introduce some functions derived from scipyin my code. Since I also need arcpy, I use the Interpreter installed inside ArcGIS, whose directory is: (C:\Python27\ArcGIS10.3\python.exe).

My issue is that, if I try to install scipy from pycharm, I get this error:

numpy.distutils.system_info.NotFoundError: no lapack/blas resources found

and they suggest me to try again to install using the prompt command. But if I execute:

pip install scipy

I get the error:

'pip' is not recognized as an internal or external command,
operable program or batch file.

I am not sure if the issue is related to the fact my interpreter is the ArcGIS one, but if yes, how can I solve the problem? Otherwise, what may the reason be?


texts: Is it possible to gray out a Google Maps map except for some area?

texts:
What i am looking for is a way to emphasize a certain area on a Google Maps map by graying out the rest of the world.

I heard about some solution of using 2 polygons of which 1 is visible and the other not and that it only works if the map is static (no map movement). However i want the user to be able to zoom and drag the map.

I am also interested in other ways i haven't thought of to emphasize an area.

update:

When i try the solution posted by radek (i.e. creating a KML with a Polygon with an outerBoundaryIs (bigger than view) and an innerBoundaryIs (being area to emphasize)) I get the problem that when i zoom in, rectangular areas are not covered by the surrounding polygon anymore. See:

Anyone with experience/solution to that?

update2:

I finally had time to test radek's solution using Google Fusion Tables. At first i had the same issues with missing tiles as with KML (see picture above), but after another load it seems like it had picked up and it was showing the tiles totally correct. I have a feeling that the jquery hack in his solution was meant for tiles that weren't loading:

setTimeout(function(){
        $("img[src*='googleapis']").each(function() {
            $(this).attr("src",$(this).attr("src")+"&"+(new Date()).getTime());
        });
    }, 5000);

As the tiles seemed to load after a refresh i haven't used the hack myself. If the bug happens for every new user I will experiement with it.


texts: Writing if/then Statements for Python Parser of ArcGIS Field Calculator

texts:
For some reason I cannot get this code to work in the Field Calculator of ArcGIS Desktop 10.6.1. I'm trying to use the Python parser to populate the field "SHAP_CLAS" with multiple conditional statements based on "AREA_LA" and "AREA_PERI_SQR".

The code I've tried so far is:

def calc(AREA_LA, AREA_PERI_SQR):
   if (AREA_LA =< 100 and AREA_PERI_SQR =< 0.65):
       return 2
   elif (AREA_LA =< 100 and AREA_PERI_SQR > 0.65):
       return 1
   elif (100 < AREA_LA < 500 and AREA_PERI_SQR =< 0.45):
       return 2
   elif (100 < AREA_LA < 500 and AREA_PERI_SQR > 0.45):
       return 1
   else:
       return 3

...and...

calc(!SHAP_CLAS!)

...in the box below the code block.


texts:
I have a problem with qgis 2.8.1 (Wien). Whenever I change the project SRC to WGS84-World Mercator (EPSG 3395), the area is calculated as zero.

How can I solve this problem? Is this a known bug?

The answer from How can I get my Qgis area calculations in metres rather than degrees? My shapefiles are projected using WGS 1984 UTM fuseau 35 Sud don't solve the problem, since the problem only appears with a metric CRS.


texts: Use an image as a fill for a polygon in GeoServer (like overlay google maps)?

texts:
Locally, I have a geoserver installed. I need to fill the polygon with a picture that is pre-cropped, along the outline of the polygon.
For example, if you use styles, and fill in one color, this is how it looks.

<PolygonSymbolizer>
        <Fill>
          <CssParameter name = "fill"> <ogc: PropertyName> #ddd </ ogc: PropertyName>
          </CssParameter>
        </Fill>
        <Stroke>
          <CssParameter name = "stroke"> <ogc: PropertyName> # ff0 </ ogc: PropertyName> </CssParameter>
          <CssParameter name = "stroke-width"> 0.5 </CssParameter>
        </Stroke>
      </PolygonSymbolizer>

It works, but I need the casting to be done with a picture, scaled, and clearly tied to the boundaries of the polygon. In short, I need the geoserver exactly as done in google maps - an example of https://developers.google.com/maps/documentation/javascript/examples/overlay-simple


texts: How to use QGIS field calculator in python?

texts:
I would like to use QGIS Field Calculator in Python to update the existing column based on certain conditions?

How can I do this?


texts: Attribute editing in QGIS on geoJson file

texts:
We have vector files in the format of geoJson files. We need to keep it in geoJson format as it is picked up by mobile devices via Restful services.

The problem is we want to change the attributes of the features in QGIS but it seems like this is not possible. We don't have to add or remove columns. We simply need to edit the existing data.

Is there a work around for this?

Are there any plugins that I can download that will allow me to do this in QGIS?


texts:
Is it possible to create a 3D pdf from qgis 3d mapview (version 3 "Girona" and above)?

I have tried with Adobe Acrobat 9 Pro Extended (3D capture) and it does not work.
The adobe works fine with other apps/programs that I use but with qgis I get an error message: "An error occured during 3D capture. Capture settings for this application may be unadapted."


texts: Keeping only 7 digit number within field name using ArcGIS Pro field calculator and regex

texts:
I have a field titled GEOGRAPHY, which contains hundreds of records containing various strings of all similar format, for example:

"Witless Bay (1001559) T 00000 (  6.5%)"
"Laniel (2485905) NO 00909 (  6.7%)"
"Contrecoeur (2459035) V 00000 (  4.9%)"

However, in order for me to join this table with my boundary layer, I need the records in this field to only contain the 7 digit number in the parenthesis, i.e. 1001559, for the Witless Bay record.

I cannot figure out how to do this through field calculator using regex. I have looked at ESRI's examples, and they only show examples for substitutions. Using LEFT() or RIGHT() doesn't work either, as the string lengths are variable.

I have tried something to the extent of the following without success with attached error:

Expression Type: Python 3

Expression:
GEOGRAPHY = update_name(!GEOGRAPHY!)

Code Block:
import re

def update_name(geo_name):
    return re.search(r"""\((.*?)\)""", geo_name)

I do not understand why this wouldn't work. The expression is valid \((.*?)\). It gives me the following error:

 The field is not nullable. [GEOGRAPHY]
 Failed to execute (CalculateField).

I understand that it gives that error because it can't return a NULL value, but should this not be returning the 7 digit number?


texts:
I'd like to use the result of an OpenStreetMap request in Overpass as a layer in QGIS project.
It would be like a VIEW on a PostgreSQL database.


texts:
I have a qgis plugin (c++) that adds a map tool inherits from QgsMapTool. My tool does nothing, juste I can activate my tool by triggering.
I test my plugin in Qgis 2.4; when I activate my tool the I passed to another tool Qgis crash and write a minimum dump file.

Dump File: qgis-20140915-111427-2760-1564-8fdd08a.dmp : C:\Users\Admin\AppData\Local\Temp\qgis-20140915-111427-2760-1564-8fdd08a.dmp
Last Write Time: 15/09/2014 11:14:28
Process Name: qgis-bin.exe : C:\Program Files (x86)\QGIS Chugiak\bin\qgis-bin.exe
Process Architecture: x86
Exception Code: 0xC0000005
Exception Information: The thread tried to read from or write to a virtual address for which it does not have the appropriate access.
Heap Information: Not Present

System Information

OS Version: 6.1.7601
CLR Version(s):

Modules

Module Name Module Path Module Version

qgis-bin.exe C:\Program Files (x86)\QGIS Chugiak\bin\qgis-bin.exe 0.0.0.0
ntdll.dll C:\Windows\SysWOW64\ntdll.dll 6.1.7601.18247
....

this is the C++ code:

MapToolAdvancedEditPolygons::MapToolAdvancedEditPolygons( QgisInterface *ptInterface )
    : QgsMapTool( ptInterface->mapCanvas () )
    , m_ptInterface ( ptInterface )
{
}

MapToolAdvancedEditPolygons::~MapToolAdvancedEditPolygons()
{
}

I initialised my plugin as below:

void GeoAdvancedEditTool::initGui ()
{
 m_ptAdvancedEditAction = new QAction ( QIcon ( ":/geoadvancededittool/images/AdvancedEditAction.png" ), tr ( "advanced edit polygons tool" ), this);
 m_ptAdvancedEditAction->setCheckable ( true );
 m_ptAdvancedEditAction->setActionGroup ( m_ptInterface->actionMoveFeature()->actionGroup() );
 m_ptInterface->digitizeToolBar()->addAction ( m_ptAdvancedEditAction );

 connect( m_ptAdvancedEditAction, SIGNAL( triggered() ), this, SLOT( advancedEditActivated() ) );

 m_ptMapToolAdvancedEditPolygons = new MapToolAdvancedEditPolygons( m_ptInterface );
 m_ptMapToolAdvancedEditPolygons->setObjectName( "GEOAdvancedEditPolygonsTool" );
 m_ptMapToolAdvancedEditPolygons->setAction( m_ptAdvancedEditAction );
 m_ptAdvancedEditAction->setEnabled(true);
}

void GeoAdvancedEditTool::advancedEditActivated()
{
 m_ptInterface->mapCanvas()->setMapTool( m_ptMapToolAdvancedEditPolygons );
 m_ptInterface->mapCanvas()->setCursor( QCursor( Qt::CrossCursor ) );
}

I don't understand why Qgis crash!


texts: Calculating the area of raster within polygon

texts: I have a raster layer with values of energy demand in each cell. I would like to know the area of the raster within each polygon (see screenshot). Any suggestions on a simple approach for this?


texts: Joining polygons to polylines based off specific point in polyline being within polygon

texts:
I am trying to join a polyline feature class and a polygon feature class. The polylines will always have at least their first point in a polygon, and each polyline consists of exactly three points. There is a chance of the last point of the polyline being in another polygon, but I do not want to join that polygon. There will be multiple polylines from within each polygon and I want to join these.

I was thinking of something like this:

arcpy.SpatialJoin_analysis(lines, polygons, linesAndPolygons, #, #, #, "POLYGON_CONTAINS_FIRST_POINT_OF_POLYLINE")

There is of course no explicit match option like above, but that is what I need. Any ideas? I was considering making a new feature class of points from the polylines, joining those to the polygons, and then joining the polylines to the points, but it seems like there should be a simpler way.


texts: How to create own OpenStreetMap tiles offline on Windows?

texts: I want to create own OpenStreetMap tiles offline on Windows. I've downloaded czech-republic-latest.osm.pbf file at geofabrik.de site, however I am unable to load it to the MapSurfer.NET and Maperative due to /16 GB/ memory overflow.


texts: Join Field Error (ERROR_000499 Table is not editable)

texts:
I'm trying to join a field from a table to a census tract shapefile based on the CT name field. I'm getting this error however:

Has anyone encountered this error before and possibly know where it's coming from?


texts: Angular - How to visualize GeoJSON on Google Maps

texts:
I am using @agm/core for visualization GeoJSON on Google Maps.

But I can't display my buildings on Google Maps.
I think the problem is the geometry that is represented as multiple аrrays.

I used FME to validate GeoJSON, and it's ok.

My data in browser:

  class:

    export interface Feature {

  geometry: Geometry;

}

export interface Geometry {

  type: String;
  coordinates: number;

}

export class Building {

  features: Feature;
  id: number;
  nameOfBuilding: string;
  nameOfCity: string;
  streetName: string;
  legalStatus: string;
  nameOfUsage: string;
}

  service:

import { Injectable } from '@angular/core';
import { Observable } from 'rxjs';
import { HttpClient } from '@angular/common/http';
import { Building } from '../classes/building';
import { map } from 'rxjs/operators';

@Injectable()
export class buildingApiService {

  constructor(private httpClient: HttpClient) { }

  getBuildings(): Observable<any> {
    return this.httpClient.get("http://localhost:8080/api/buildings")
  };
}

  app.component.ts

import { Component, OnInit, Input } from '@angular/core';
import { buildingApiService } from './services/building.api.service';
import { Building } from './classes/building';
import { extend } from 'webdriver-js-extender';
import { AgmCoreModule, AgmDataLayer } from '@agm/core';
import { reduce } from 'rxjs/operators';

@Component({
  selector: 'app-root',
  templateUrl: './app.component.html',
  styleUrls: ['./app.component.css']
})
export class AppComponent implements OnInit {

  lat: number = 45.2461012;
  lng: number = 19.8495081;
  zoom: number = 15;

  constructor(private _buildingApiService: buildingApiService) { }

  buildings: any = [];

  ngOnInit() {

    this._buildingApiService.getBuildings()
      .subscribe
      (
        data => {
          this.buildings = data;
          console.log(data);
          })
  }

}

  app.component.html

<agm-map [latitude]="lat" [longitude]="lng" [zoom]="zoom">

  <agm-data-layer [geoJson]="buildings">

  </agm-data-layer>

</agm-map>

How to solve it?


texts: Why do I get unexpected coordinates after converting EPSG 102014 to EPSG 4326?

texts:
I have this Europe map (EPSG 102014) and I am trying to convert it to EPSG 4326 using geoconverter.hsr.ch

After the conversion is done, I expect that coordinates are in the range of -180,180 and -90,90, but instead, I get this range: 10.00,10.10 and 30.00,30.10.

If you check this example here, you can see that the countries have the same shape, it's all looking great, but they are all concentrated in a tiny area of Africa.

How can I convert between the two EPSGs and get the right latitude and longitude?


texts: How can i make drainage density map by using drainage map

texts: I made drainage map, but I don't know how can make drainage density map, I just heard that I have to convert the lines into points and then use HeatMap, but I did and my results do not make sense. Who can help me to figure out it?


texts: Get QgsLayoutItemMap from existing QPT Layout

texts:
I created a QPT model with a map inside QGIS. The problem is when I try to select the map under pyQGIS :

qpt = QDomDocument()
qpt.setContent(content)
layout = QgsLayout(QgsProject().instance())
layout.loadFromTemplate(qpt, QgsReadWriteContext())
mapItem = layout.itemById('Map') # 'Map' is the name of the item map

mapCanvas = iface.mapCanvas()
mapItem.zoomToExtent(mapCanvas.extent())

It returns this error :

AttributeError: 'QgsLayoutItem' object has no attribute 'zoomToExtent'

I know there is a class QgsLayoutItemMap but I don't know how to get it from the QgsLayoutItem map.

And I have another question : in QGIS composer, there is a button when I select my Map called "Set to map canvas extent". What is the equivalent in python ? mapItem.zoomToExtent(mapCanvas.extent()) and then mapItem.refresh() ?


texts: Importing google maps/bing points: which projection?

texts:
I am trying to use data imported through one of the google maps/Bing API, and am confused about the projection. Using the APIs, data comes under a non-geographic format, such as json or xml, try:

http://maps.googleapis.com/maps/api/geocode/json?address=Davis&sensor=false

It seems that the projection of google maps is epsg:3857. However, according to this stack post and this page, things are more complicated:

  All of this further confused by that fact that often even though the
  map is in Web Mercator(EPSG: 3857), the actual coordinates used are in
  lat-long (EPSG: 4326).

So what should I do, when I import a point, and want it, say in the final projection WGS 84 (epsg:4326)? I tried:

Declaring as 3857, then converting to 4326: gives me a point out of space (see point_3857_to_4326 in code below: -0.00109, 0.00034)
Declaring directly as 4326: seems to make sense (see point_4326)

Where is the issue? That the points themselves are not in 3857, or that they are, but conversion is not well done? This post indicates in chapter "the solution" how to do the conversion, but I just don't get it...

Example with R:

# not run: library(ggmap)
# geocode("Davis", output="latlon")

library(sf)
point <- st_point(c(-121.7405167, 38.5449065))

point_4326 <- st_sfc(point, crs=4326)
point_3857_to_4326 <- st_sfc(point, crs=3857) %>%
  st_transform(4326)

point_4326
point_3857_to_4326


texts: How to add time and elevation dimension attributes to Geoserver rasters?

texts:
In geoserver, when I edit a postgis raster layer and click on the Dimensions Tab, I can't enable either of the Time or Elevation dimensions. The following message appears next to each of the respective enable checkboxes:

Time

  Cannot enable, no attribute of type Date found

Elevation

  Cannot enable, no attribute of type Number found

Where do I add these attributes? I do have a date column in my raster postgis table.


texts: Determining projection of Texas county shapefile?

texts:
I have a Shapefile of San Jacinto County in Texas.  It thinks it is NAD83 UTM Zone 15N, but the dataset isn't jiving with any other dataset.

Max X: 1,073,721   Max Y: 11,223,970
Min X: 895,728      Min Y: 11,013,769

I don't usually work with datasets from this part of the world, so I don't know which projections are in common use.  Anyone who is familiar with the area know what projection these extents might fit?

This is the county I'm talking about:


texts: GDAL create Grid brings invalid Grid on the first run

texts:
I am playing around with gdal and numpy to create a random grid in QGIS. The following code sniplet creates an grid an tries to load into QGIS. But the the result layer is invalid. If I run the script a second time, right after the first run, the layer is created and loaded without any problem...do I miss something in my code sniplet?

from osgeo import gdal
import numpy as np

fn = "d:/newgrid.tif"
rasterband = np.random.random((10,10))
driver = gdal.GetDriverByName("GTiff")
ds = driver.Create(fn, xsize=10, ysize=10, bands=1, eType=gdal.GDT_Float32)
ds.GetRasterBand(1).WriteArray(rasterband)
rlayer = iface.addRasterLayer(fn)

The error message on the first run is:  Invalid Layer : GDAL provider Cannot open GDAL dataset d:/newgrid.tif:          `d:/newgrid.tif' not recognized as a supported file format. Raster layer Provider is not valid (provider: gdal, URI: d:/newgrid.tif


texts: Overlapping line aggregation+trim in PostGIS, while summing values

texts:
I'm sorry for the vague title, I didn't really know how to name this. I am working on a project that visualizes and simulates 17th century trade. For this project we also use so called Minard lines to show how certain goods are distributed. To clarify I made this screenshot:

However, as you can see every route is a seperate line (Amsterdam <> Jakarta and Rotterdam <> Jakarta for instance) and they all overlap with each other. This does give a nice Tron effect, but we actually want to merge them. We want to merge the lines that overlap (while adding their values), and segmentize those that don't (while remaining their values).

I've been thinking a lot about this on how to do this in postgres and postgis, but I'm having real troubles with it. Functions like ST_union aggregates everything, Could anyone point me in the right direction? Here's a screenshot of the database:


texts: SRTM DEM data appear mismatching in BW/Saturation or even brightness when gathered

texts:
I'm trying to make a DEM data of Algeria as I tried to download about 15 SRTM data but when added on ArcGIS and gathered they appear with diffrent saturation/brightness each, I tried then again to join the rasters with the Mosaic to raster tool, it actually matches but decreases in information detail and harshly in brightness..
https://ibb.co/QCwb50z


texts: How to return a polygon geometry from overlapping rows that share a same value in another column

texts:
This one is not very easy to explain but basically I have a table with postcodes and a table with area boundaries.

These boundaries have an area_id and postcodes that match the postcode table

This area_id repeats itself for each postcode that the boundary contains, so the geometries should overlap at times, for example:

Area| postcode

A | NW1

A | NW2

A | NW3

B | NW2

B | SE4

Illustrated, something that would look like the below, with the middle polygon overlapping since it belongs to both area A and B:

Query used so far and doesn't work is

create table test as
select area_id, st_union(b.geom) as geom
from postcodes b, area_boundary c
where b.postcode = c.pcd
group by c.area_id


texts: Calculating sum of values in point shapefile within moving window using QGIS

texts:
I have a point shapefile that has a 'rank' attribute with numbers between 0-100. I need to calculate the sum of this value for points within a 10km circular moving window anywhere within my area of interest. The result would be a raster, where each pixel represents the sum value of the ranks taken from all point features within a 10km radius. In ArcGIS Desktop, I could do this using the Point Statistics tool.
Is there a way of doing this in QGIS?


texts:
I operate behind a firewall at work, and so to access the plugins and other external features for QGIS I need to use a proxy.
We also have our own WMS server - which doesn't like the proxy.
To access the internal WMS, I need the proxy off, to access an external WFS I need it on.

Is there a way to control the proxy (as found in Settings > Options > Network)?


texts: Finding out what geographies ACS 2015 5 year variables are reported at for use with Census API?

texts:
I'm trying to query the US Census 2015 5 year ACS API for certain variables at the block group level, but unfortunately not all variables/tables are reported at the block group level.

How can I find out what geographies a particular variable or table is reported at?

I've looked all over and can't find any documentation on this.


texts:
I am using a .dll file methods in mapbasic

declare method openmultiplefilesdlg Class "mbextensions",lib "mbextensiondlgs.dll"
   (declaration variable) as string

but it shows error:

  Found [method] while searching for [sub]

and for function declaration

declare function splitline(byval input as string, byval unit_id as integer) as integer

it shows error:

  Expecting sub or function


texts: How to create graduated Pie Charts from CSV in QGIS 2.0.1 Defour?

texts: I am trying to make a map in QGIS 2.0.1 ( running on my Mac OS X Maveriks) with the proportionally sized (graduated) pie-charts. I have started with the csv conversion to kml just to place my sample sites on the map, and I am wondering if I can just modify this csv adding some extra columns indicating i) sample size and ii) percentage of each sector for the pie chart. Or there is much less complicated way?


texts: Folium/Leaflet display marker without radius?

texts:
I would like to plot a point in a Folium map. Not a circle, just a point (i.e. without a radius).

I would like for them not to zoom in if I do, but just stay as a point.

I have

   fmap.circle_marker(location=[row["lat"],row["lon"]], popup=row["tower"])

I've tried setting the radius to 0, and the fill_color/fill_opacity to 0 as well, but it still displays a circle.


texts:
I'm trying to replicate the following example from package getCRUCLdata:

# Create a raster stack of temperature from tmp
# files in the tempdir() directory.
download.file(
   url = "https://crudata.uea.ac.uk/cru/data/hrg/tmc/grid_10min_tmp.dat.gz",
   destfile = file.path(tempdir(), "grid_10min_tmp.dat.gz")
)
CRU_tmp <- create_CRU_stack(tmp = TRUE, dsn = tempdir())
CRU_tmpcode here

However, all that I get is an error message: 'gzip' not recognized as an internal or external command, operable program or batch file.
Can anybody tell me how to fix this?

I'm running on:

sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17134)

Matrix products: default

locale:
[1] LC_COLLATE=Portuguese_Brazil.1252  LC_CTYPE=Portuguese_Brazil.1252
LC_MONETARY=Portuguese_Brazil.1252 LC_NUMERIC=C
LC_TIME=Portuguese_Brazil.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] R.utils_2.9.0      R.oo_1.22.0        R.methodsS3_1.7.1  getCRUCLdata_0.3.0
raster_2.9-5       sp_1.3-1

loaded via a namespace (and not attached):
[1] Rcpp_1.0.1        lattice_0.20-38   codetools_0.2-16  digest_0.6.19
hoardr_0.5.2      rappdirs_0.3.1    grid_3.6.0        R6_2.4.0
data.table_1.12.2 rgdal_1.4-3       tools_3.6.0
[12] compiler_3.6.0

Alternatively, I would like to know how to convert the .dat file, obtained by decompressing the file downloadable from the link in the example, to a raster stack.


texts: GIS Clustering with drawing predefined square area around points as cluster boundaries

texts:
I have a layer of points and I would like to create clusters where a predefined square area of 5km by 5km is the boundary for the cluster.

Is there a way to iteratively draw square polygons around the points so that the polygons are not overlapping and contain at least 1 point and not more than a 100 points? The polygons can be rotated to get the best fit.

It would be better not to start with one cluster containing the maximum number of points as a starting point but rather to try evenly distribute the points in the clusters especially with dens areas with lots of points. The cluster size of 1 should be only created if there is no possibility to combine at least two points.

I am imagining the following rules.
start with:
- find the biggest cluster of points where area is bigger than 5km by 5km, draw the boundaries around it so it can be in next step divided into even 5km by 5km squares(big one can be rotated to get the best fit so minimizing the empty space)
- find the next cluster by size and repeat the procedure
- find any number of points which can fit the 5km by 5km
- in the last instance find the single points and draw 5km by 5km squares

It would be like drawing a customized grid.

I had to create this task very quickly so I used square buffers and join by location algorithms but had to adjust dens locations by hand. I will have more points (around 20k) in the future and would like to automate the procedure.

I used QGIS but I was thinking this could be probably done also with R.

How could this be achieved and what tools could be used to do it?


texts: Formatting CSV file to export it as DXF?


texts: Mutli feature polyline into individual features polyline

texts: Is there a tool which can split mutli feature polyline into an individual feature polyline, without spliting polylines that cross other polylines in the same shapefile?


texts: Creating multi-ring buffers (donuts) on a feature collection GEE

texts:
I have a large feature collection that I am creating a multi-ring buffer on using GEE. My problem is that the .difference algorithm only works on a feature and not a feature collection. I tried to flatten the feature collection but that is still not working. How can I do this?
See code snippet below
////A function to add buffer rings to the refugee sites, recommended buffer size 40-50km
var ROI_a=RefugSites.map(function(ft){return ft.buffer(10000)});
var ROI_b=RefugSites.map(function(ft){return ft.buffer(30000)}).flatten();
//trying to use the inner circle here to subtract it from the second buffer
var ROI_b=ROI_b.map(function (f){return f.difference(ROI_a)});

var ROI_c=RefugSites.map(function(ft){return ft.buffer(40000)});
//and to do the same thing here, for the third ring/donuts
var ROI_c=ROI_c.map(function (f){return f.difference(ROI_b)});


texts: Error: shape mismatch when plotting certain DEM Geotiffs in 3D - using 'plot_surface' from Matplotlib

texts:
I'm trying to plot various DEM geotiff rasters in Python 3.7 using plot_surface in matplotlib  but encountered an error that I can't solve. I have 2 DEM rasters, the small one plots correctly, the other doesn't.
Here is the plot of the small raster:

Here is the error I get with the larger raster:
line 30:
surf = ax.plot_surface(X,Y,Z, rstride=1, zorder = 0, cstride=1, cmap=plt.cm.RdYlBu_r,  linewidth=0, antialiased=False)

..\axes3d.py", line 1504, in plot_surface
    X, Y, Z = np.broadcast_arrays(X, Y, Z)

..\stride_tricks.py", line 258, in broadcast_arrays
    shape = _broadcast_shape(*args)

..\stride_tricks.py", line 189, in _broadcast_shape
    b = np.broadcast(*args[:32])

ValueError: shape mismatch: objects cannot be broadcast to a single shape

The two geotiffs are downloadable here.
Here is the script:
import numpy as np
import gdal
import geopandas
from geopandas import GeoSeries
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from shapely.geometry import shape

fig, ax = plt.subplots(figsize=(8,5), subplot_kw={'projection': '3d'})
plt.ioff()

demFile = r"C:\Temp\small_DEM_8bit.tif"     # works
demFile = r"C:\Temp\DEM_8bit.tif"           # fails

ds = gdal.Open(demFile)
gt  = ds.GetGeoTransform()
dem = ds.ReadAsArray()

xres = gt[1]
yres = gt[5]
X = np.arange(gt[0], gt[0] + dem.shape[1]*xres, xres)
Y = np.arange(gt[3], gt[3] + dem.shape[0]*yres, yres)

X, Y = np.meshgrid(X, Y)
Z = dem

print("X: {}. Y: {}. Z: {}".format(len(X),len(Y),len(Z)))

surf = ax.plot_surface(X,Y,Z, rstride=1, zorder = 0, cstride=1, cmap=plt.cm.RdYlBu_r,  linewidth=0, antialiased=False)
print("plotted surface")
ax.set_zlim(0, 100)
ax.view_init(30,-105)
ax.set_zlim3d(-1.01, 200)

fig.colorbar(surf, shrink=0.4, aspect=20)
fig.show()

Both rasters are 8 bit unsigned geotiffs, both produced by gdal with the same extract process. Obviously the second raster is larger but still only 681KB and the dimensions of the X,Y,Z match each other for both the small and large tif. I tried transposing the Z (np.transpose(Z)) but it didn't work either.
Please can someone put me out of my misery as to why I get the error?


texts: Set layer in edit mode in QGIS processing

texts:
Set a layer editable through the Python console works fine here (QGIS 2.18.4):

layer=iface.activeLayer()
if not layer.isEditable():
    qgis.utils.iface.actionToggleEditing().trigger()

But I do not know how to set a input layer editable in a QGIS processing script:

##Input_layer=vector
# ...
inlayer = processing.getObject(Input_layer)
if not inlayer.isEditable():
    #???

Can anybody help me with this?


texts: Displaying x,y coordinates not working correctly in ArcGIS Desktop?

texts:
My issue is similar to Displaying X, Y coordinates on ArcGIS not working but the solution there did not work for me.

I am importing Excel data into ArcMap, and have already saved it as a CSV. When I Display X, Y Data... from the dropdown list, it shows the points at "0,0" on the map. When I change the projection from the WGS84 to Africa_Albers_Equal_Area_Conic, it shifts the points to center on the African continent.

When I brought the CSV into ArcMap, it converted my formatting of decimal degrees (correctly, I believe) and then I padded them with zeroes as that seemed to be other peoples problems (pictures).

The points themselves are all clustered within 1.0m of each other. Again incorrect, but there seems to be slight variation in x, y.


texts:
I have a set of polygons and I want to find the bbox for each polygon using sf.
With st_bbox I obtain the bbox for the whole set of features but I want to obtain bbox for each feature.

I try to use group_by() but the output is always the same.

library(sf)

g <-  st_sfc(st_point(1:2), st_point(3:4)) %>% st_buffer(dist = 1)
s <-  st_sf(a=3:4, g)
s %>% group_by(a) %>% st_bbox()%>% st_as_sfc


texts: Rasterize a SpatialPolygonsDataFrame: Error in if (x2a < rxmn)

texts:
I want to rasterize a SpatialPolygonsDataframe derived from a tessellation algorithm and get the following error:

Error in if (x2a < rxmn) { : missing value where TRUE/FALSE needed

The data for creating the tessellation comes from a random percolation map and contains only 0 and 1. The particular thing here is, that the error seems to appear only for a certain arrangement for the 0 and 1 - and it is not quite obvious to me what is causing it. It appears that .polygonsToRaster assumes that the polygon is showing holes in its geometry?

The code to reproduce the error:

library(dismo); library(deldir); library(igraph)

m <- matrix(c(1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0,
  0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1,
  0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1,
  0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0,
  1, 1, 0, 0), ncol=10, nrow=10)

r <- raster::raster(m)

# Cluster identification (clustering of adjoining pixels) ----
clusters <- raster::clump(r, direction = 4)
# Convert cells with values to points
cluster_points <- raster::rasterToPoints(clusters, spatial = TRUE)
# Create voronoi polygons
tess <- dismo::voronoi(cluster_points)

# Convert to raster
newr <- raster::raster(nrow = 10, ncol = 10, ext = raster::extent(c(0,1,0,1)))
tess_r <- raster::rasterize(tess, newr, field ='clumps', fun = "mean")

Any ideas on how to fix it or another way to rasterize the polygon?

EDIT
Spacedman solved the error, but his solution leads to another problem:

The initial number of cells is 10x10 - after his rasterization, there are 11x11 cells. When I now try to crop the spatialpolygonsdataframe to an extent of c(0,1,0,1) (the same as the raster has) and try rasterize:

randomcluster_spdf <- raster::crop(r, extent(0, 1, 0, 1))
# Convert to raster ----
randomcluster_raster <- raster::rasterize(
  tess,
  raster::raster(nrow = 10, ncol = 10, resolution = c(1/10, 1/10), ext =
  raster::extent(c(0,1,0,1))),
  field = tess@data[, 1],
  fun = "mean"
)

**Error in if (x2a < rxmn) { : missing value where TRUE/FALSE needed**

... is back. Is there the possibility to rasterize and keep the original dimensions of 10x10?


texts: Copying Network Analysis Layer and Service Area Polygon in ArcPy?

texts:
I have built a geoprocessing tool using arcpy modules in python. Screenshot shown below:

In this tool I am creating Service Area Polygons and selecting census boundaries within those Service Area Polygons to extract the demographic data. The part where I am creating the service area polygon and copying them to output is shown below:

#PART 3: CREATING SERVICE AREA POLYGON USING DRIVING DISTANCE AS COST FACTOR AND INTERSECTING WITH CENSUS LAYER AND CREATING OUTPUT

Service_Area = "Service Area"

# Process: Make Service Area Layer
arcpy.na.MakeServiceAreaLayer(road_network, Service_Area, "Length", "", distance)

# Process: Add Facility to the Network
arcpy.na.AddLocations(Service_Area, "Facilities", "facility", "", "")

# Process: Solve the Network
arcpy.na.Solve(Service_Area)

# Copy Polygons to Service Area output
arcpy.CopyFeatures_management(Service_Area, service_area_output)

# Process: Create a new Feature layer - Census layer
arcpy.MakeFeatureLayer_management(cen_blk,"blocks_lyr")

# Process: Create a new Feature layer - Census layer
arcpy.MakeFeatureLayer_management(service_area_output,"service_area_output")

# Process: Select Layer By Location...
arcpy.SelectLayerByLocation_management("blocks_lyr", "HAVE_THEIR_CENTER_IN", "service_area_output")
if census_output=="#" or census_output=="":
    census_output==None

else:
    # Process: Copy the intersected census layer to output FC and add it into current workspace
    arcpy.CopyFeatures_management("blocks_lyr", census_output)

gp.AddMessage (" ")

I cannot seem copy the service area polygon to a feature class. I am getting a following error message.

Is CopyFeature not the right method to copy the service area polygons?


texts: Georeferencing multiple images as GeoTIFF (same size and extents) in QGIS

texts:
I have 50 images and I want to georeference them all on to the exactly same place. They have all the exactly same size and extends. I've begun to load the same georeferencing points for every single image and georeferenced it... but this takes a way too long.

Is there a easy way to georeference them with one coordinate file as batch all to the same place?


texts: Set transparent raster values with Python

texts:
I have loaded a few raster layers in QGIS 2.18.

I select a few of them. Next, I run the script:

The first option

#Change properties selected layers
selectedLayers = iface.legendInterface().selectedLayers(True)
for layer in selectedLayers:
    print layer.name()
    pixel = QgsRasterTransparency().initializeTransparentPixelList(255,255,255)
    layer.renderer().setRasterTransparency(pixel)
    layer.triggerRepaint()

 print "Finish"

Property of the selected layers has not changed.

The second option

#Change properties selected layers
selectedLayers = iface.legendInterface().selectedLayers(True)
for layer in selectedLayers:
    print layer.name()
    rasterTransparency = layer.renderer().rasterTransparency()
    listPixels = []
    pixel = QgsRasterTransparency.TransparentSingleValuePixel()
    pixel.red = 255
    pixel.green = 255
    pixel.blue = 255
    pixel.percentTransparent = 100
    listPixels.append(pixel)
    rasterTransparency.setTransparentSingleValuePixelList(listPixels)

#Refresh map
for layMap in qgis.utils.iface.mapCanvas().layers():
    layMap.triggerRepaint()

print "Finish"

Property of the selected layers has not changed.

What am I doing wrong?


texts: Compute pixel area of a large raster in Google Earth Engine

texts:
I want to compute a total area of a large raster dataset using .reduceRegion(), but receive the "Earth Engine capacity exceeded" Error. Is it possible to overcome it? For example create a vector grid over an AOI, map .reduceRegion() over each sell and then sum up the results?

The code is following (also available on Earth Engine):

var aoi = ee.Geometry.Polygon(
        [[[1.0167888602549056, 32.24107677376602],
          [1.0167888602549056, 3.3274515320624296],
          [42.457218547754906, 3.3274515320624296],
          [42.457218547754906, 32.24107677376602]]], null, false);

var dataset = ee.Image('JRC/GSW1_0/GlobalSurfaceWater');
var selected_data = dataset.select('occurrence').clip(aoi);

var area = selected_data.reduceRegion({
      reducer: ee.Reducer.sum(),
      geometry: aoi,
      scale: 15,
      maxPixels:1e13
    });

print(area);

Map.addLayer(selected_data, {}, 'Target image');
Map.addLayer(aoi, {}, 'AOI');
Map.centerObject(selected_data);


texts: exclude value from feature analysis classification (ArcGIS extension)

texts:
I am using feature analysis to classify small (mm) round features from the landscape. I have calculated the circularity of each shape and I want to exclude those which are not round enough.

Does anyone know how to exclude specific values from a classification using feature analysis?

I'm also using ArcMap 10.2


texts: How to add GeoServer Vector tiles in React using OpenLayers API?

texts:
We are creating a React application and we want to display vector tile layers from GeoServer. The basemap is required from ArcGIS REST service using the following code and it works properly
const basemap = new TileLayer({
      source: new XYZ({
        attributions: 'Tiles © ArcGIS',
        url: 'https://server.arcgisonline.com/ArcGIS/rest/services/' +
            'Canvas/World_Dark_Gray_Base/MapServer/tile/{z}/{y}/{x}'
      })
    });

As for the vector tile layers, we use Mapbox Vector Tiles (MVT) and the code is shown as the following but it did not work.
    const initialStateSource = new VectorTileSource({
      url: 'https://smartcommunityhealth.ahc.umn.edu/gwc/service/tms/1.0.0/SOLAP%3Aus_states_simplified@EPSG%3A900913@pbf/{z}/{x}/{y}.pbf',
      format: new MVT(),
    });

    const initialStateLayer = new VectorTileLayer({
      visible: statesVisible,
      source: initialStateSource,
      style: geoOutline
    });

How can I solve this problem?


texts: Optional Parameter Not Being Optional in a Toolbox

texts:
I have been creating a toolbox that contains 4 parameters. Two of my parameters involve editing the labels of a shapefile. I want to have these two parameters as optional: if the user does not want to have the labels edited, then don't edit them!

I should mention that I'm not up to that point yet though. The two parameters are the label text size and a boolean to see if the user wants the text to be bold. I have set the parameters to be optional in my getParameterInfo method:

def getParameterInfo(self):
    """Define parameter definitions"""

    # First parameter = layer with customized features
    param0 = arcpy.Parameter(
        displayName = "Source layer with customized features",
        name = "source",
        datatype = "Feature Layer",
        parameterType = "Required",
        direction = "Input")

    # Second parameter = layers to which customized features will be applied
    param1 = arcpy.Parameter(
        displayName = "Layers with features to be updated",
        name = "updlayers",
        datatype = "GPValueTable",
        parameterType = "Required",
        direction = "Input")

    param2 = arcpy.Parameter(
        displayName = "Label Font Size",
        name = "labelFontSize",
        datatype = "GPDouble",
        parameterType = "Optional",
        direction = "Input")

    param3 = arcpy.Parameter(
        displayName = "Bold Font:",
        name = "labelBoldFont",
        datatype = "GPBoolean",
        parameterType = "Optional",
        direction = "Input")

    param1.columns = [['Feature Layer', 'Features']]

    params = [param0, param1, param2, param3]

    return params

However, although they are set as optional (param2 and param3), the tool still requires an input in both cases. Here is a screenshot of what I see when I open my toolbox. Notice the green bubble next to the last two parameters, which indicates that they are required:

Is there a way for the last two parameters to be truly optional? As in, the user does not have to do anything with them if they don't want to? I'd imagine that this is something I can do in the getParameterInfo() method, but I am unsure of how to do this. Would I be able to set default values for these parameters in some way?

One extra thing to mention is that if I click on the check box for bold font to enable it, and then click it again to disable it, the tool will run. This is a temporary workaround. However, it is not ideal. I would hope to have the user not have to interact with these options if they don't want to.


texts: Creating composite raster from rasters with different spatial extent without losing information using ArcGIS Pro

texts:
I am creating a composite raster using several individual rasters with different spatial extents (captured from different times), see image below.

The idea here is that I want to create one mosaic composite raster without losing any information (i.e. calculating means of repeating raster values) in order to have a more complete picture. I would normally use the Mosaic to New Raster tool for a situation like this, however, in this particular case there is quite a bit of overlap between the rasters so I just want to make sure I don't lose any data.


texts: Speedup function to calculate distance from centroid of building to river segment in PostgreSQL

texts:
I am trying to measure the distance from the centroid of a building to a river segment using the following function.

CREATE OR REPLACE FUNCTION river_foo_mapper () RETURNS trigger LANGUAGE plpgsql AS $$
BEGIN
new.river_distance := subquery.distance

                    FROM (with center as (
                        select st_transform(st_centroid(new.geometry),3857) as geom from
                        foo )
                        select st_distance(a.geom,st_transform(b.geometry,3857)) as distance
                        from center as a, rivers as b
                            order by st_transform(b.geometry,3857) <-> a.geom
                            limit 1) AS subquery ;
  RETURN NEW;
  END
 $$;

If I run the SQL as a select statement it takes less than 2 seconds but if I run it as a function it takes forever to finish. How best can I optimize this and is my function correct.


texts:
I have a layer of pipes,water conductors, that I merged  under the condition of having the same material and they're touching each other, it did a good work but  some features which  have the same material but not touching each other are also merged.
the code below shows the query I executed to merge the features :

drop table if exists touches_material;
CREATE TABLE touches_material AS
SELECT ROW_NUMBER() OVER() AS id,
       sub_query.*
FROM (
    SELECT array_agg(DISTINCT a.id) AS old_id, a.material,
           ST_Union(a.geom) AS fusion
    FROM pipe AS a,
         pipe AS b
    WHERE a.id <> b.id
      AND a.material = b.material
      AND ST_touches(a.geom, b.geom)
    GROUP BY a.material
) AS sub_query

So I was thinking about finding a function that clears the merge when the touch is not realized and the pipes have the same material. How can I achieve this goal?


texts:
I am using QGIS 3.10.
I created raster map indicating the cost of traversing each cell in the north-south and east-west directions which I have used to create a cumulative cost map with the module r.cost. It works fine.
However, I try to use the same friccion map with the cost allocation tool and it doesn't work, I get error, although I use the same files than with the cumulative cost tool.
Anyone know if the parameters to use cost allocation map are different? Any suggestions?


texts: What software is available on OSX for working with DEMs?

texts: I have several adjacent DEM files and would like to combine them into a mosaic and then export to DEM and TXT. I know a Windows solution is to use MICRODEM but before installing Windows 7 I would like to know if there is a Mac (OSX) based option.


texts:
I am running a QGIS network script called "Transit-tools-QGIS" on QGIS 2.12.3

I am interested in this script for my network:

https://github.com/chourmo/QGIS-Transit-tools/blob/master/scripts/Shortest%20paths%20-%20car.py

I get the error:

No module named VectorWriter.

Can anyone advise me on this issue? I am also a qgis newb.


texts: Multigeometry relations to single geometries, through Overpass-Turbo to KML

texts:
I'm trying to extract relations from Openstreetmap, using the Overpass-Turbo service, saving as a KML-file to be viewed on a suitable mobile phone App. I need the resulting KML to have simple geometries for it to work properly on the mobile app.

Relations in OSM can have multiple geometries (e.g. more than one polyline within one relation). An export of data, extracted from OSM using the Overpass-Turbo service, respects this. So the raw KML will be structured like this:

<Placemark>
      <name>Some name</name>
 <MultiGeometry>
        <LineString>
          <coordinates>1234567890...</coordinates>
        </LineString>
        <LineString>
          <coordinates>

Proper GIS software should be able to cope with this. But KML-files are increasingly used for viewing on mobile Apps, which may expect the kml-file to contain only single polylines.

The app 'Maps.me' for example, displays such multi-polylines without gaps between what should be separate parts of the polyline (often making a real mess).

What would be the simplest way to extract data like this from OSM using Overpass, if we need to convert such multi-geometries into simple polylines (or polygons etc) so that they display on less sophisticated apps correctly?

Obviously one solution is to bring data into GIS software, using this to manipulate the data, but I'd prefer a solution which avoids the use of GIS software completely.

I was aiming to distribute an Overpass-Turbo query to a small number of non-GIS using colleagues, enabling them to use this to extract specific data from OSM themselves to view on their own mobile phones.


texts: How to create a plain boundary from a large parcels shapefile

texts: Basically I just want the outline of a shapefile, and nothing inside it.  I have tried creating a buffer and dissolving but i still get some stuff in the interior of the polygon.  I just want the outside.  Also tried clipping the parcel file with itself but that leaves all the parcels intact.  There are a lot of edges on the outside so creating a trace polygon is not gonna work for me.  Is there another way?


texts: How to create a label combining different font sizes or types?

texts:
I would like to place the elevation-number of a point in a different front size and centred under its name:

Is that possible?

(That is my actual labeling: label || '\n' || elevation)


texts: Permanently Re-arrange fields of attribute table (feature class and tables)?

texts: I am trying to Find out simple Geoprocessing Tool that can physically re-arrange the fields of attribute table (for both feature class +  tables), I have found some posts here, but nether of them have provide a good and a simple tool that can do the  job.


texts: How to remove QGIS polygon boundaries when using graduated style?

texts:
When I use a graduated style and I set the outline of the polygons to transparent I can still see the border. If I set the width to "Hairline" and/or the style to "No pen" I can still see them.

I have found multiple similar questions with very manual intensive solutions that are simply not practical in my case:

Removing border of symbol in QGIS?
Forcing QGIS to remove polygon borders?
QGIS renders outlines where none should be

This seems to be an old feature request, am I missing something?

I'm using QGIS 2.18.


texts: What Does the Cell Value in a Composite Raster Mean?

texts:
When I make a composite raster in QGIS using the GRASS r.composite tool and then click on a certain cell I get these values in the blue/green/red/composite rasters:
blue = 124
green = 124
red = 172
composite = 15861.

I know that the blue/green/red represent the brightness values of the given bands. What does the 15861 represent and how is it calculated?


texts: AttributeError: MapDocObject: Unable to save on mxd.save/saveACopy?

texts:
I'm trying to use arcpy to save an .mxd and am having an issue.  It doesn't matter if I use mxd.save or mxd.saveACopy I get the same error.  I know it doesn't look like I am doing anything w/ it here, but it's just a portion of my code:

import arcpy
import arcpy.mapping as map
arcpy.env.overwriteOutput = True

env.workspace = r"C:\MGIS\FinalProject2.gdb
mxd = map.MapDocument(r"C:\MGIS\Finalproject2.mxd")
df = map.ListDataFrames(mxd)[0]

for df in map.ListDataFrames(mxd):
df.scale = 2000000

#Add layers into map document
addLayer = map.Layer(OKXcounties)
addLayer1 = map.Layer(contourLabelText)
addLayer2 = map.Layer(contoursClipLyr)
addLayer3 = map.Layer(clipOKXLyr)
map.AddLayer(df, addLayer, "AUTO_ARRANGE")
map.AddLayer(df, addLayer1, "AUTO_ARRANGE")
map.AddLayer(df, addLayer2, "AUTO_ARRANGE")
map.AddLayer(df, addLayer3, "AUTO_ARRANGE")

# Add contour Labels
for lyr in map.ListLayers(mxd):
    if lyr.name == "contourLabelText":
        lyr.showLabels = True
    if lyr.name == "Contour Features":
        lyr.transparency = 100

mxd.saveACopy(mxd)

I get the following error:

Traceback (most recent call last):
  File "C:\Python26\ArcGIS10.0\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 312, in RunScript
    exec codeObject in __main__.__dict__
  File "C:\MGIS\geog485\FinalProject2\mxdtest.py", line 7, in <module>
    mxd.saveACopy(mxd)
  File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\utils.py", line 181, in fn_
    return fn(*args, **kw)
  File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\_mapping.py", line 668, in saveACopy
    self._arc_object.saveACopy(file_name)
AttributeError: MapDocObject: Unable to save.  Check to make sure you have write access to the specified file and that there is enough space on the storage device to hold your document.

Space is not the issue, I have 60 GB of free space and I've checked the windows properties of the .mxd file.  Read-only box is unchecked and I've "allowed" all permissions for all groups:

Does anyone have any idea how if this is a bug or how I can fix this?

I'm using ArcInfo 10.0 SP4.


texts: Extract information of a raster overlayed on another raster with a lower resolution directly in python

texts: I have two aligned rasters with the same extent, one has a resolution of 5 meters and the other of 125 meters. I would like to place the 5 meter one over the 125 meter one and assign to these smaller pixels the value of the large pixel on which they land on. How can I do this in pure python code? The raster I/O is not a problem.


texts: How to force caching of image layers in Leaftlet?

texts:
I'm working on a weather map project with Leaflet. I have several weather parameters as PNG overlays. For each parameter there's 60 PNG files, representing different points in time. I have implemented a simple routine which will allow the user to step forward and backward through time like this:

function getUrl(parameter,time) {
    return '/charts/leaflet/'+parameter+'-'+time+'.png';
}

function setUrls(time) {
    console.log(time);
    levels.setUrl(getUrl('levels',time));
    capesfc.setUrl(getUrl('capesfc',time));
}

var increaseTime = L.Toolbar2.Action.extend({
    options: {
        toolbarIcon: {
            html: '<i class="fa fa-step-forward" />',
            tooltip: 'increase time'
        }
    },
    addHooks: function () {
        if (timeIndex < 160) timeIndex++;
        else timeIndex = 100;
        setUrls(timeIndex);
    }
});
var decreaseTime = L.Toolbar2.Action.extend({
    options: {
        toolbarIcon: {
            html: '<i class="fa fa-step-backward" />',
            tooltip: 'decrease time'
        }
    },
    addHooks: function () {
        if (timeIndex > 100) timeIndex--;
        else timeIndex = 160;
        setUrls(timeIndex);
   }
});

As you can see, my system is based on directly changing the layer URL which will trigger Leaflet to reload the layer from the new URL. Unfortunately, as the files are about 400k each, that creates somewhat of a delay and flicker. I was planning to include some animation function and/or a quick scroll through some slider and for that it's to stuttery.
Is there a way to force the caching of the images or a smarter way to handle things alltogether, e.g. when a layer is selected, Leaflet could preload all layers or something like that?


texts: Extracting coordinates from line Shapefile

texts:
I am a complete loss when attempting to use QGIS to produce a csv file of points as latitude and longitude from shapefile data so that I can plot them on gores for a small world globe. All my attempts have resulted in just a list of consecutive numbers in a column. The data I have chosen has 250 points and that will be ideal for my purpose, a 55mm diameter globe to illustrate a talk on Mercator and the Ordnance Survey

I am using the data from http://openstreetmapdata.com/data/generalized-coastlines The shapefile set is the set from coastlines_z1 ... with extensions .dbf .prj .shp .shx I have tried to follow the instructions for QGIS in this reference,

1) Add a new layer
2) use the file coastlinesz1.shp
3) toggle editing
4) open the layer’s attribute table

At this point the table has only one column with the numbers 0 to 249 in it.

Somewhere along this process I have obviously made at least one mistake, please can anyone tell me where?.

I am sure that more will be needed, but there is a plethora of information out there and deciding what is relevant is beyond me at the moment.


texts: Open feature class in in_memory workspace ArcGIS Pro SDK

texts:
I am aware that one can open a feature class in a geodatabase using the code below.

Geodatabase geodatabase = new Geodatabase(new DatabaseConnectionFile(new Uri("path\\to\\sde\\file\\sdefile.sde")))
FeatureClass table = geodatabase.OpenDataset<FeatureClass>("FeatureClassName");

But how can I open a file that's saved in memory such as "in_memory/testFile"?

I have combed Esri's resources to no avail.


texts:
If I get a map back from a WMS server is the scale always constant over them map regardless of the reference system?

So if I have a map which is 100 pixels wide, will the point at 50 pixels always be half way between the minimum and maximum values of the sides of the bounding box that defines the map?

So can I always linearly interpolate the map coordinates based on the pixel position?

Or can there be cases where 50% of the way across the map in pixels is not 50% of the way between the minimum and maximum of the width of the bounding box for some reason, like the projection is not constant?


texts: Atlas Generation output from QGIS is returning empty files

texts:
I have set up a project and print composer to generate an atlas so I can print off some 300 pages of feature coverages.

I have selected "Preview Atlas" and entered in settings for Atlas Generation and Controlled by Atlas from the item properties.  I am able to navigate through the pages as normal, but when I "Save as Image", QGIS returns the empty files with the names of the coverage layer; which doesn't work.

Have I missed any steps or settings? I also downloaded and installed QGIS 2.14.3 Essen, which behaves the same.


texts: Can geoserver handle normal get request?

texts:
I am relatively  new to the geoserver, backend - frontend GIS architecture. I know how to generate a WFS with geoserver from my PostGIS database and retrieve this data in my openlayers frontend.

However, I still would like to know how I can handle the case of a mobile app, sending attributes via an URL with gps coordinates etc. to the postGIS database through Geoserver. The app generating the URL already exists. I just would like to know if Geoserver can handle a basic get request from this URL to send the data to the database. This would happen every 30 sec or so.

EDIT:
It is not a WFS which is generated by the app. so i have to create the PostGIS geometry within Geoserver with the available coordinate values from the URL


texts: QGIS for generic surface documentation setting a custom user scale and CRS

texts:
I would like to use QGIS for the documentation of generic flat surfaces (so I don't need coordinates relative to a sphere but to a cartesian flat plane) but just overlays some images of surfaces with known dimensions and adding information to them. I would like to open a raster file and set the origin to 0 and the corresponding points to a chosen reference systems so that the units of measurement are right and I can perform some estimation of area and length in the new reference system.

I'm trying to set a user-defined projection system using proj4. I'm looking for a suitable CRS for a 2D surface.
At the moment I am opening the raster using the georeferenced and adding some controls to place it on the qgis project. Is there any better approach for using QGIS for documenting generic flat surfaces?


texts: How to convert a "Polygon" to a "GeometryCollection" type in PostGIS

texts:
With this simplified example (In reality I have thousands of points):

CREATE TABLE test
(
  geog geometry(GeometryCollection)
);

INSERT INTO test VALUES (ST_ConcaveHull('MULTIPOINT(-109 29, -110 28, 30 10, -109 29)',0.99));

I get the following error:

ERROR: Geometry type (Polygon) does not match column type (GeometryCollection)

How can I convert the Polygon type to a GeometryCollection type?

ST_Multi and ST_Collect only convert to MultiPolygon , not GeometryCollection.


texts: Adapting a ConvexHull fonction to a ConcaveHull one

texts:
I have a database with over 150 000 points in a EPSG 2950 projection. As I wanted to cluster those points in polygons, I used a ConvexHull fonction. The goal was to create clusters where I would find AT LEAST 5 points at a MAX DISTANCE of 75 meters from one to another.

The fonction I used was the following:

Select *
INTO new_table
 FROM (Select ST_ConvexHull
    (unnest
     (ST_ClusterWithin(geom, 75)))
      AS geom from my_table) sq
Where ST_NPoints(geom) >= 5

The fonction worked really well but now I would like to try to have clusters where the perimeter is a bit more representative of the points inside (The convexHull can leave big gaps where there is no points).

To do so I tried replacing the ConvexHull fonction for a ConcaveHull one:

Select *
INTO new_table
 FROM (Select ST_ConcaveHull
    (ST_Collect(unnest(ST_ClusterWithin(geom, 75), 0.70)))
      AS geom my_table) sq
Where ST_NPoints(geom) >= 5

This specific function returned the following error...

  ERROR:  the function unnest(geometry[], numeric) doesn't exist
  LINE 4:     (ST_Collect(unnest(ST_ClusterWithin(geom, 75), 0.70)))

I have tried some variants (as using the ConcaveHull function before the FROM clause) but the closest I got was getting 1 ConcaveHull over ALL my points and therefore ignoring the constraint of having AT LEAST 5 points at a MAX DISTANCE of 75 m.

Any ideas to solve my problem?


texts: Displaying QGIS layer after loading with Python

texts:
I'm successfully loading several layers with Python, but can only get to see them in the map window if I manually zoom to extent afterwards.

Once the last layer is loaded I call zoom to extent in the code on the layer, but it seems not to have an effect. Ending the loading with iface.mapCanvas().refresh() similarly have no effect.

I've looked at these two related Q&As:

Load and display a layer from SQL Server using PyQGIS and
https://gis.stackexchange.com/questions/80910

but not found an answer to why I can't see the layer. What am I missing?


texts: How to analyze and combine overlaping elevation maps into one?

texts:
I have 30 overlapping maps that contain tree density heights (measures between 1-100%) of a large area, and I need to get the highest point by pixel and create a new map with all these highest points if it makes sense.

Basically, say in point (pixel) "A" in map 1, 2...30 the maximum elevation is 25 m, and in the next pixel "B" (in all the other maps) is 20. I would need to create a new map (a new matrix) that condenses all these elevations/densities (and will look like 25, 20, x, y, etc).

Someone has suggested making a program to analyze all of the maps by pixel, but I'm new to this, so any help, or suggestion will be appreciated.

Ok, based on a comment I got, I'll add more detail here about the data:
I have the data in three different formats (the same data) these being tif, .dtm. and ASCII. I have a section of a river, and the 30 maps overlap on top of each other exactly. They have measurements of three density (0, where there is no "vegetation" in the pixel, to 100, all is vegetated), now each of these overlapping maps contain this same density (0-100) for different heights, namely 1-30m.
For example: in map 1, I would have calculated already from another program the points between the "bare ground" and 1m, so my map "1" contains everything ABOVE 1m, my map 2, will contain everything after 2m, etc until reaching 30m.

From these maps, I'll need to produce a map with the highest three height by pixel.

Thanks


texts: Grouping values after Unique symbology classification in ArcGIS for Desktop?

texts:
I have a layer that I am using the unique values, many fields symbology category for but I have two different values that I would like to have displayed with the same symbol and label for that symbol in the legend.

If I didn't have a legend I would just assign them both the same symbol and not worry about it but with the legend it would end up having the same thing twice on it that way.

Is there a way to make a value in this situation be some thing like -1,* where as long as the first value field is -1 it doesn't matter what the second one is?


texts: Loading GeoWebCache tiles into OpenLayers

texts:
I have geotiff image that I want to show on a site locally (when there is no internet connection and browsing to localhost).

I've seen that OpenLayers 3 has XYZ layer but the format doesn't match the GeoWebCache output folder.

How can I view tiles locally from GeoWebCache... or is there free app which can generate XYZ tiles from Tiff ?


texts: QGIS editFormConfig setNotNull not working in Python

texts:
When trying to use editFormConfig(), I get setNotNull is undefined error shown below.

   AttributeError: 'QgsEditFormConfig' object has no attribute 'setNotNull'

Below is the part of the code used.

  layer.editFormConfig().setNotNull(idx, True)

Is setNotNull in python binding? How could I fix the issue?


texts: Automatic search of most northern/eastern/western/southern polygons in sets of polygons

texts: I have a large set of round polygons, each of them is divided into 4 equal quarters by 2 perpendicular lines with azimuths of 45 and 135 degrees. Is there a way to proceed with an automatic selection of most northern/eastern/western/southern polygon in each of these sets of four?


texts:
I'm using GeoNetwork 3.0 where users can upload their own layers from the interface to my GeoServer. GeoNetwork 3.0 uses OpenLayers3. For my application, I need to determine the crs from this layer so I can transform it to EPSG:3857, if it's not already.

How do I get the projection from a layer in OL3? Or can I tell GeoServer to transform all uploaded layers to EPSG:3857 by default?


texts: How can you link a qtdesigner form from the QGIS itself, and to call a code made from the console?


texts: Selecting raster values at random locations in a raster in python/arcgis

texts:
I have a large dataset of integer rasters, and was hoping to find a way to generate a "random" selection of pixels for four of the rasters, while making sure all 5 classes, which I've specified based off the pixel values, are equally represented with the selection. Each raster's values ranges from 0-200, with 0-100 representing percent tree coverage and 200 representing water (aka data I want to ignore for this case). I already added a new field to the raster called "TreeClass," and filled that field with numbers 1-6, where 1-5 are different levels/classes of tree coverage and 6 is water/no data values that I want to ignore. In all for each raster, I would like 10 pixels, 2 in each class 1-5.

In addition to randomly selecting these pixels, I would like to record their locations (plus values and classes, if possible), so that I can do a time series analysis on the selected pixels using the other rasters

Any suggestions on the best method for this? I have thought about several approaches:

Converting the four rasters to polygon based on TreeClass and then use the Create Random Points tool with the polygon as a constraint. However, after using the Raster to Polygon tool, it doesn't appear to have come out right, and is very slow. Additionally, I don't see a way in which this tool would allow me to store the pixel location (not the coordinates, but samples, lines location data of the pixel).
Converting the raster to numpy array and then using python tools to generate random pixels (not yet sure how I would go about this). That would make it easy to store the x, y location of the pixel and then use them to get the same pixel location's value for the other numpy arrays (arrays that I would convert from raster). However, I don't see a way to convert the raster to a structured array so that not only the pixel value would be present and accessible but the class value too. Perhaps I could convert to numpy array that stores the pixel value, somehow add a dimension to the numpy array, and then calculate the class value in python to store into the new dimension.

I'm hoping there's a simple solution and I'm just over-complicating things, as I tend to do. But I wanted to get some feedback that will help lead me in the right direction.

Thanks so much in advance.


texts: Quantifying amount of sunlight falling in a region

texts:
I am trying to objectively quantify or provide some form of indication that could point toward the amount of sunlight

Anyone knows if there are any datasets (either instrumental or remotely sensed data) or derivational approaches that provides such information on the amount of sunlight received?


texts: QGIS Server working well on local VM but remote server config fail

texts:
Trying to understand how to setup QGIS Server I started on a local Ubuntu VM. After a day of debuging I got something that's work.

I try to replicate the same on my hosted VM (GCP). I got the classical  :

<ServiceExceptionReport xmlns="http://www.opengis.net/ogc" version="1.3.0">
  <ServiceException code="Service configuration error">Service unknown or unsupported</ServiceException>
</ServiceExceptionReport>

When I hit the server IP + cgi-bin/qgis_mapserv.fcgi. So QGIS Server itself work.

But when I try to get the capabilities I got a 500 error, which lead me to server issue more than QGIS Server.

/cgi-bin/qgismapserv.fcgi?SERVICE=WMS&REQUEST=GetCapabilitie

Here the VirtualHost that I have.

<VirtualHost *:80>
  ServerAdmin my@email
  ServerName qgis.test # -> I have binded the ip and the name to my /etc/hosts the time I have a real URL to add to my DNS.

  DocumentRoot /var/www/html

  # Apache logs (different than QGIS Server log)
  ErrorLog ${APACHE_LOG_DIR}/qgis.demo.error.log
  CustomLog ${APACHE_LOG_DIR}/qgis.demo.access.log combined

  # Longer timeout for WPS... default = 40
  FcgidIOTimeout 120

  FcgidInitialEnv LC_ALL "en_US.UTF-8"
  FcgidInitialEnv PYTHONIOENCODING UTF-8
  FcgidInitialEnv LANG "en_US.UTF-8"

  # QGIS log (different from apache logs) see https://docs.qgis.org/testing/en/docs/user_manual/working_with_ogc/ogc_server_support.html#qgis-server-logging
  FcgidInitialEnv QGIS_SERVER_LOG_FILE /var/log/qgis/qgisserver.log
  FcgidInitialEnv QGIS_SERVER_LOG_LEVEL 0

  FcgidInitialEnv QGIS_DEBUG 1

  # default QGIS project
  SetEnv QGIS_PROJECT_FILE /home/sysadm/01_General_test.qgs

  # QGIS_AUTH_DB_DIR_PATH must lead to a directory writeable by the Server's FCGI process user
  FcgidInitialEnv QGIS_AUTH_DB_DIR_PATH "/home/sysadm/qgisserverdb/"
  FcgidInitialEnv QGIS_AUTH_PASSWORD_FILE "/home/sysadm/qgisserverdb/qgis-auth.db"

  # See https://docs.qgis.org/testing/en/docs/user_manual/working_with_vector/supported_data.html#pg-service-file
  SetEnv PGSERVICEFILE /home/sysadm/.pg_service.conf
  FcgidInitialEnv PGPASSFILE "/home/sysadm/.pgpass"

  # Tell QGIS Server instances to use a specific display number
  FcgidInitialEnv DISPLAY ":99"

  # if qgis-server is installed from packages in debian based distros this is usually /usr/lib/cgi-bin/
  # run "locate qgis_mapserv.fcgi" if you don't know where qgis_mapserv.fcgi is
  ScriptAlias /cgi-bin/ /usr/lib/cgi-bin/
  <Directory "/usr/lib/cgi-bin/">
    AllowOverride None
    Options +ExecCGI -MultiViews -SymLinksIfOwnerMatch
    Order allow,deny
    Allow from all
    Require all granted
  </Directory>

 <IfModule mod_fcgid.c>
   FcgidMaxRequestLen 26214400
   FcgidConnectTimeout 60
 </IfModule>

</VirtualHost>

Any idea ?


texts: Can a DXF File contain raster images (TIFF/JPG)?

texts:
I am maintaining a large Program with a graphical editor. Our Program can export DXF and it can  show Background Raster Images of various types (TIFF, JPG, BMP etc.).

Is it possible to embed these Raster Images in the DXF export so that they will be displayed when such a file is opened with AutoCad?


texts: Create QGIS 3.14 PostGIS local connection without user name or password?

texts:
I used to be able to create local PostGIS connections in QGIS (as of 3.12.x) without supplying the user name or password. On the local test machine, PostgreSQL is set up to trust authentication, so it will accept any credential local (and the default user name -- the login user).
However, this seems to break down after upgrading to the new QGIS 3.14 LTR a few days ago.
When trying to create a default local connection, I was always asked to provide username and password (and service name). None of these were required previously.
Is there a way to work around this and get the old behavior back?
(This is with Ubuntu 20.04 LTS)


texts: How do I show floors to on a polygon using CesiumJs

texts: I have a gltf 3D building model on cesium. I would like to color code the range of floors. i.e. 2-5 floor, red 5-10 blue etc. How do I go about doing this?


texts: Set-up GRASS for use in QGIS 3.10 with Mac OS X

texts:
I want to use GRASS 7.6 with QGIS 3.10 installed on my Mac (10.13 - High Sierra).

I have installed QGIS 3.10 for Mac via this link : https://qgis.org/en/site/forusers/download.html
I have installed GRASS 7.6 for Mac via this link : http://grassmac.wikidot.com/downloads

QGIS and GRASS both work separately, however when I am unable to use GRASS via QGIS. After activating the plugin in QGIS, I do not have access to any of the features (New Mapset, Open Mapset, etc.). Accessing at Plugins > GRASS > GRASS Options > General, I get the following error message :

Currently selected GRASS installation is not valid

I have the application located in my Applications directory.

What are the steps that I should do to be able to use GRASS via QGIS?


texts: I have recently begun accessing Planet data through the Planet Explorer interface.  I am most interested in the Planet Basic imagery.  These are the images that have not yet been orthorectified.  For nearly all scenes I see the Analytic and Visual products, but not the Basic.  Is there something I need to do to make this option available?


texts:
I followed this post, the first answer, to clip a USGS Quad in order to remove the borders so that I can tile them. Basically it utilized the lat/long at the corners in a .txt file to import a rectangle using "add delimited text layer". I then reproject the rectangle to match the CRS of the USGS quad. After I clip it I get black edges, see image.

I offered the error as a comment in the aforementioned post and it was answered that its likely from re-projecting the rectangle and that I should use -dstalpha or -dstnodata to make it white. When I use those two options, I just get a colored rectangle. I actually want that area transparent for when I tile multiple USGS Quad maps. I even tried to re-project the quad map to the CRS of the rectangle.

Does anyone know how I can clip the USGS quad without having the black corners for tiling pruposes? FYI... I am using QGIS 2.2 on Windows 7.

EDIT:

The mask that I have compiled using the delimited text layer matches the area that I would like to clip perfectly. So I dont think that I have an error. I downloaded the Newburgh and  the Cornwall on Hudson Quads. Reading the corner coordinates I developed the following two .txt files:

Newburgh

  Nr;WKT
  1;POLYGON((-74.00 41.5, -74.125 41.5, -74.125 41.625, -74.00 41.625, -74.00 41.5))

Cornwall on Hudson

  Nr;WKT
  1;POLYGON((-74.00 41.375, -74.125 41.375, -74.125 41.5, -74.00 41.5, -74.00 41.375))


texts: Can't merge Polygons - Error: "layers must all be of the same geometry"

texts:
I have several polygons (all of which have the same CRS and attributes) and I want to merge these into a single layer.  I had hoped this would be a quick, easy process in QGIS but the merge functions (listed under data management tools and the one that is part of MMQIS) fails and instead report that the layers must be of the same geometry.

I am aware that Grass may provide a solution to this issue - as per  Error in merging 2 shape files of different geometry types - but my knowledge of QGIS is limited and I would rather focus on that than needing to get my head around GRASS as well. Is there a simple function in QGIS that can resolve this issue?


texts: Get coordinate value displayed at the buttom of qgis window with pyqgis

texts: Does anyone know how I can get the coordinate value displayed at the buttom of the qgis window? I have been searching through the qgis api but am not able to find anything suitable..


texts: How to compare Coronavirus deaths between countries

texts:
Leaving aside the problem of how deaths are counted (in some countries a doctor diagnoses the case, others insist on a positive test, some only count hospital deaths), how do we accurately compare between countries?

Raw numbers of deaths cannot be used because some countries have larger populations than others. Deaths per 100,000 population is a little better but in the UK, the majority of the deaths are in the large cities of London and Birmingham. These large death figures concentrated into two cities are diluted by the uninfected population elsewhere. This website https://www.bbc.co.uk/news/52311014 suggests comparing like with like e.g. London with New York, both being globalized cities, but how do you define "like"?

My idea is to calculate number of deaths per 100,000 in each region in the country (e.g. Public Health Area, Superoutput area - as small a region as data can be obtained for). Then create a raster of death/100k per grid square in each region. International comparisons could then be made by ranking countries by number of grid squares in each category. That is those with large numbers of high value death/100k squares would be considered most impacted by the virus. Deaths would be thus compared to the immediately surrounding population rather than the population of the entire country.

Is this a reasonable approach?


texts:
I'm a new on GIS, And I'm coding a program that can interact on WMS map, I using the https://github.com/ClemensFischer/XAML-Map-Control library.

I have two WMS capabilities links.

https://www.wms.nrw.de/geobasis/wms_nw_dop?request=GetCapabilities
https://geoportal.regionalverband-braunschweig.de/geoserver/rrop08/wms?request=GetCapabilities&version=1.3.0&service=wms

With the first WMS capabilities link, I can load the image right way with config:

WMS URL: https://www.wms.nrw.de/geobasis/wms_nw_dop?
WMS Layer: WMS_NW_DOP
CRS: EPSG:25832

The Map Can Show Smooth.

But with the second WMS Capabilities, the config I using is

WMS URL: https://geoportal.regionalverband-braunschweig.de/geoserver/rrop08/ows?
WMS Layer: RROP08 Verkehr Luft
CRS: EPSG:25832

the image is totally white.

And I found some different between two Capabilities XML. With the second link. The name of layer have space. and it is longer than the first.

And I think that, seem I not understand about the WMS capabilities XML file, so how to read the XML file by the right way?

What is real WMS URL?, What WMS layer should be selected? and more more thing about the XML tag's mean?


texts: Point count within borders, in two different layers, using QGIS

texts: I have a layer of georeferenced points, and a layer of vectorized borders. Is there a way I can use QGIS to calculate how many points fall within each area of the borders layer?


texts:
Possible Duplicate:
  Automated correction of parcel geometry based on area measure

I need a function to modify area polygon in PostgreSQL.
The major problem is automated correction of parcel geometry based
on area measure, so I need to modify area polygon in PostgreSQL


texts: Speeding up Leaflet with GeoJSON data app?

texts:
I am working on a Leaflet app which will be displaying five layers with about 1k points each.

For our purpose (basically in order to locate most of the points near to a road) we don't want to use clustering; instead, we want to filter the amount of points on each zoom level, so there are never too many points on the map at any given time (say no more than a few hundred).

I have tried two different approaches on this issue. First, filtering the points on each zoom level according to a field on the GeoJSON data (Z_level).

Second has been creating different GeoJSON files, one for every possible zoom level we'll be using.

I've found that filtering is a little bit (10-20%) faster than layer reloading, but am finding a big hit in performance when zooming in after zoom level 13 or 14 when there are a few thousands of points loaded (not rendered). Even more surprising for me, the biggest performance hit is when zooming out from level 14 to 13 (about 8 seconds) and even more when zooming from 13 to 12 (almost 12 seconds!!).

The rest of zoom transitions are taking some 2 to 3 seconds, which is acceptable.

Here is the basis of the code I'm using for each case:

1) Filtering:

    // Filtering points for each Layer

var towns = L.geoJson(null, {

   filter: function(feature) {
        if(feature.properties.Z_Level <= mapZoom)
            return true;
      },
  pointToLayer: function (feature, latlng) {
    return L.circleMarker(latlng, {
      title: feature.properties.Name
       })
  }
});

map.on("zoomend", function(e) {
            map.removeLayer(towns);
            towns = L.geoJson(data, {
               filter: function(feature) {
                    if(feature.properties.Z_Level <= mapZoom)
                        return true;
                  },
               pointToLayer: function (feature, latlng) {
                    return L.circleMarker(latlng, {
                        title: feature.properties.Name
                        });
                    }
                }

            map.addLayer(towns);
            });

And here is the second one (layers):

// Different layers for each Zoom Level

for(i=minZoom; i<=maxZoom;i++){

  towns[i] = L.geoJson(null, {

  pointToLayer: function (feature, latlng) {
    return L.circleMarker(latlng, {
      title: feature.properties.Name
       })
  }
});
}

map.on("zoomstart", function(e) {
    mapOldZoom = map.getZoom();
 });

map.on("zoomend", function(e) {

    mapZoom = map.getZoom();

       if (map.hasLayer(towns)) {

                map.removeLayer(towns[mapOldZoom]);
                map.addLayer(towns[mapZoom]);
                }
 });

Can anyone send any ideas on which would be the best approach? Am I losing other ways to improve speed on a case like this?