texts: how to calculate the ndvi % greater than 0.5 in a particular zone and also cloud percentage

texts:
i have a problem whem i am working to calculate the ndvi from landsat sattellite of rainy season i have assumed the pixel value above 15000 are clouds and calculate NDVI.then i am assigning the clouds pixel in ndvi 1.

The main question is how to calculate zonal  min max value of ndvi and percentage of cloud and percentage of ndvi greater than 0.5 in a particular zone.(zone are the vilage boundary of District).
please help to solve this problem in arc gis

i have already calculated mean max and min value of ndvi but i want to calculate polygone(zone) wise NDVI >.5(%) and also same in clouds.


texts: R: discrepencies between pointDistance and costDistance

texts:
I'm trying to do some landscape analyses using R, and am having issues where using the pointDistance function from the raster package and the costDistance function from the gDistance package can give significantly different results depending on the raster properties. Here is a trivial example:

library(raster)
library(gdistance)

x <- c(5, 370 , 900)
y <- c(5, 370 , 900)

coords <- cbind(x,y)

pd <- pointDistance(coords, lonlat = FALSE)

#works; dist_matrix matches pd
r <- raster(ncol=1000,nrow=1000, crs=NA, xmn=0, xmx=1000, ymn=0, ymx=1000)
r <- setValues(r,1)

tr <- transition(r, mean, directions=8)
tr1 <- geoCorrection(tr, type="c", multpl=FALSE)

dist = costDistance(tr1,coords)
dist_matrix = as.matrix(dist)

#doesn't work
r <- raster(ncol=100,nrow=100, crs=NA, xmn=0, xmx=1000, ymn=0, ymx=1000)
r <- setValues(r,1)

tr <- transition(r, mean, directions=8)
tr1 <- geoCorrection(tr, type="c", multpl=FALSE)

dist = costDistance(tr1,coords)
dist_matrix = as.matrix(dist)

When the raster has spatial resolution of 1 for each pixel, it works, but anything else, it starts to fail. In my case, it's resulting in a 1-10% difference between each site, which is a lot. I'm adding barriers (lakes) to my rasters which have less of an impact...

So how do I deal with this? Depending on the extent I'm working with, my rasters range from tens of megabytes to hundreds of megabytes in size, with each pixel representing a 50x50 meter area. Eventually, when I have the entire process ironed out, it's going to drop to a 10x10 meter area for each pixel, which is going to be the limit of what I can work with computationally. I won't be able to get it down to a 1x1 resolution like the trivial example that works.

Edit: Turns out that the trivial example only works when the x component and y component of a single point are the same value. Also, changing spatial resolutions of the raster doesn't seem to affect the costDistance results, so at least that function seems consistent regardless of the size of the pixels in the raster.


texts: Is it possible to set the permissions Geoserver creates new directories with?

texts:
I'm trying to automate the process of pushing from a server to Github as a way of tracking Geoserver changes, and pulling from Github as a way of deploying new work. The push side works fine, but if I create a new layer in Geoserver, and then try to edit that elsewhere and pull changes back to the server it fails.

Steps:

Publish new layer in geoserver
Push that from the server, pull it to my own machine
Make any changes - e.g. adding height.ftl or description.ftl
Push that from my machine, try to pull it back on the server.

At this point, I get a "permission denied" error, and the issue seems to be that Geoserver isn't giving these directories group write permissions; specifically it's setting: drwxr-xr-x

A manual chmod g+w fixes this, but I'm trying to get the process more automated.  Is there a way to set the permissions with which Geoserver creates these directories?

[adding with later edit: it occurred to me to check the umask settings, and that seems to be set to 0002 for all of my user account, the user Geoserver runs as, and the user Github runs as.  They're also all members of each others' groups. If I've understood this correctly, it means the server's default behaviour is to give group write permission to new directories and files, so Geoserver must be explicitly withholding that permission.]


texts:
I am getting strange vertical lines on my hillshade, anyone experienced anything similar?

I am very new to ArcGIS.


texts: QGIS layer blending: keep polygon outline un-blended?

texts:
Version 3.8.3.  QGIS' layer blending is a spectacular feature that I use often.  For example, when displaying a polygon layer over a satellite image background I'll set the polygon layer blend mode = "multiply".

However, this blending also affects the polygon boundary.  I'd like the boundary to remain un-blended (i.e. blending mode = "normal").  Is there a way to do this without displaying a second, hollow copy of the polygon layer on top with its blending set to "normal"?

In summary, how do I simultaneously set polygon fill = blend multiply, with polygon stroke = blend normal?


texts:
I am trying to convert a SHP vector file of X,Y,Z bathymetry. I am trying to use the GDAL Rasterise (vector to raster) tool.

When I run it I get the following error in the log. I have tried looking for the /bin/sh: gdal_rasterize: file on my Mac but can't locate it.

QGIS version: 3.8.0-Zanzibar
Qt version: 5.11.2
GDAL version: 2.4.1
GEOS version: 3.7.1-CAPI-1.11.1 27a5e771
PROJ version: Rel. 5.2.0, September 15th, 2018
Processing algorithm…
Algorithm 'Rasterize (vector to raster)' starting…
Input parameters:
{ 'BURN' : 0, 'DATA_TYPE' : 5, 'EXTENT' : '473840.0,479842.0,6564270.0,6570650.0 [EPSG:23031]', 'FIELD' : 'field_3', 'HEIGHT' : 2.5, 'INIT' : None, 'INPUT' : '/Users/willhutchins/Documents-Offline/2mbathy.gpkg|layername=2mbathy', 'INVERT' : False, 'NODATA' : 0, 'OPTIONS' : '', 'OUTPUT' : '/Users/willhutchins/Documents-Offline/2m bathy.tif', 'UNITS' : 1, 'WIDTH' : 2.5 }

GDAL command:
gdal_rasterize -l 2mbathy -a field_3 -tr 2.5 2.5 -a_nodata 0.0 -te 473840.0 6564270.0 479842.0 6570650.0 -ot Float32 -of GTiff /Users/willhutchins/Documents-Offline/2mbathy.gpkg "/Users/willhutchins/Documents-Offline/2m bathy.tif"
GDAL command output:
/bin/sh: gdal_rasterize: command not found

Execution completed in 0.10 seconds
Results:
{'OUTPUT': '/Users/willhutchins/Documents-Offline/2m bathy.tif'}

Loading resulting layers
The following layers were not correctly generated.

/Users/willhutchins/Documents-Offline/2m bathy.tif

You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.


texts: Why dissolve creates many extra features?

texts: I have a shp file in Arcmap 10 that has ~16600 polygons. Some overlap. I want any polygons that are overlapping to be dissolved into single polygons. I want non-overlapping polygons to remain the same. When I use Dissolve as suggested here How to dissolve overlapping polygons? I get about 93000 polygons. Why is it touching my non-overlapping polygons? Only thing I can guess is that it's something to do with polygons that were multipart polygons to begin with and now I've unchecked "create multipart polygons." If that's it, how can I work around that? If not, what gives?


texts: what is the Maximum Number of replica download possible from download.geofabrik.de?

texts:
I am trying to update more than 5 months of Nominatim europe data from download.geofabrik.de. For 8 to 10 days its works fine but after that it showing error
failed to open stream: HTTP request failed! HTTP/1.1 500 Internal Server Error

Is there any restriction in downloading replication from download.geofabrik.de server per IP?


texts: Classifying fields and getting statistics in ArcMap?

texts:
I am using ArcMap 10.6.

I have a table with country data. I am interested in Asia only, so I created a new layer (and shapefile). Several countries have missing data, which was filled with the same value (-50). One of the fields is Life Expectancy and another is Access to Drinking Water(population percent). I classified the Life Expectancy field (about 30 unique values) into 5 classes (in the symbology tab). I don't know how to group the countries by class (in the attribute table) so I can calculate various stats by class.

For example, What is the mean Access to Drinking Water, by Life Expectancy class?


texts: How to analyze census data and display on map (shapefile) in ArcGIS

texts:
I have shapefiles representing a country, states of that country or towns etc... For each shapefile there is associated data with them.

The country shapefile has 100 rows and maps to .csv files with 100 rows, same goes for towns, states etc. Basically it gives me a flexibility at how much level I can experiment with.

Now, let's say I chose the states data and I wanted to find whether folks with education are more likely to get jobs than those who don't, the first and easy part was to join the shapefile to its .csv file that contains data about employment.

After that I simply got confused.. the CSV files contains raw data and how to use them confuses me, an example (fake data):

SID STYPE   SNAME        MALEP  FEMALEP EMPLOYED  UMEMPLOYED EDUCATION
S1  STA     Alabama      4345   1003    100000    2000       1000
S2  STA     Alaska       1456   8343    8888      4000
S3  STA     Arizona      90000  1247    1000
S4  STA     Arkansas     18831  1213    0

How do I use the data above? I don't think a select by attributes will do me any good as it shows just a blue border, and I don't know how I could possible query that? It doesn't give me any specific data, it is just an aggregate.

I think my best option is to use either categories or quantities symbology.

However, I'm not sure how to get meaningful data. If I was to use quantities (I think quantities is the way to go), do I simply normalize people who are employed against people having education?
How would I display the output.. density dots or gradient symbology/color?


texts: Help binding Dynamic Map Layer popup in Esri Leaflet?

texts:
I've been experimenting with a web map using Esri Leaflet examples and I've run into a bit of a problem getting my popups to load. I'm able to successfully load the data from the REST service and it displays properly, however no popup is created when I click each feature.

Sites = L.esri.dynamicMapLayer("http://gis.montgomerycountymd.gov/ArcGIS/rest/services/General/Fiber_PLOI/MapServer",{});
Sites.bindPopup(function (error, featureCollection) {
if(error || featureCollection.features.length === 0) {
  return false;
} else {
  return 'Name: ' + featureCollection.features.properties.NAME;}
});


texts: Calculating buffer distance for different groups of parallel lines in QGIS

texts:
I have a number of parallel lines (see here for the background: How to extend lines to Bounding Box in QGIS?). My final aim is to create the Centerline for each group of lines as illustrated here:

What I tried: I want to group these lines together in such a way that each line is grouped together with it's closest neighboring line. Then I want to create a dissolved buffer around these groups. See the next screenshot to better understand what my intention is: in this case, the lines should be grouped together to four groups. The left- and rightmost of the lines of each group should become the boundary of a polygon.
So the question is: how can I process the single parallel lines in way to get polygons around each group of lines?

Shapefile is available on the following link: https://drive.google.com/drive/folders/1FI7cPVd4qKUU3d8C1_vNj_RfWQ1HqZqX?usp=sharing


texts: Loading Sentinel 2 images downloaded with QGIS in SNAP software?

texts:
I have downloaded Sentinel 2 images using the semi-automatic classification plugin in QGIS, but I cannot open these files in SNAP software.

I can load .jp2 products in SNAP but I obtain a error when I try to open the XML file in SNAP.

  Error: No appropriate product reader found

My SNAP version is the 5.0 and the image identifier is: S2A_OPER_MSI_L1C_TL_SGS__20161117T180146_A007341_T30TXN


texts: Map location of unknown point with distances to two known points

texts:
I'm trying to develop a function in R that will locate a point (A) on a map given distances between that point A and two other points (B and C), for which exact coordinates are known. The problem I'm coming up against is that my approach is sensitive to the shape of the triangle, and thus the number of conditions in my algorithm is higher than I think it needs to be, probably because my middle-school maths aren't good enough, so I'm looking for somebody to be really condescending, tell me how silly I'm being, and help me to a general approach to this problem.

For example, I'm looking for a single algorithm that can place point A in all of these circumstances:

(x,y): B(1,2); C(4,1), distances: A to B = 6, A to C = 2
(x,y): B(1,1); C(4,2), distances: A to B = 2, A to C = 3
(x,y): B(1,1); C(4,1), distances: A to B = 2, A to C = 6

For the record, I'm aware that there are two possible points in space that satisfy these assumptions. I'm happy to have both.

Lastly, I'll point out that while this sounds like homework, it's actually ecological research, albeit performed by somebody apparently underqualified to be doing it.


texts: Creating raster layer from dictionary in PyQGIS?


texts: Select features using a selected feature in ArcGIS Pro

texts:
Using a selected feature is it possible to select from a separate file to then extract?

In this case I'm trying to extract points from within a given parish using a whole country parish layer.


texts: Batch exporting MXD to PDF files using ArcGIS Desktop?

texts:
Nico Burgerhart made a fantastic script called "Batch export MXD to PDF" back in 2008.
    http://arcscripts.esri.com/details.asp?dbid=14872
    Any ideas how to implement it in ArcGIS 10? Something using Python perhaps?
    Data Driven Pages does not have anything comparable that i could find.

Here's the steps in read me text file which doesn't jive with ArcGIS 10's new menu bar:

Tool: Batch export MXD to PDF
Purpose: Saves all MXDs in the selected directory to PDFs in the selected output folder
Author: Nico Burgerhart (nicoburgerhart@hotmail.com)
Date: 31 Jan. 2007
INSTALLATION NOTES
------------------
1. Open ArcMap
2. Select Tools > Macro's > Visual Basic Editor
3. Select File > Import file
4. Import BatchExportMXDToPDF.bas
5. Select File > Close and Return to ArcMap
6. Select Tools > Macro's > Macro's
7. Select the BatchExportMXDToPDF mactro
8. Click Run


texts:
I have just created a date field attribute.

In Proprieties > Fields of the layer I made this:

And the attribute tables shows this:

What I want to see is the date in this format: 07-01-2019 and not 2019-01-27.

What do I need to do to have a permanent solution?


texts: I'm trying to export a shapefile to a dwg-file to use it in AutoCAD. When I export the shapefile using the "Export to CAD" tool, the dwg-file that's being created is not georeferenced (but I need it to be). Anyone who can help me?


texts: Geoprocessing tool not executing when Map's SpatialReference is Web Mercartor

texts:
I am facing a strange issue when trying to execute GP service when my map's spatial-reference is Web Mercator. I appreciate all help.

I have a JavaScript API application. Initially I only had my map service used in this application. The Map Service I am using is in NAD27. I execute a GP service from JavaScript API and it works fine. In fiddler I can see the request being send to the server.  Now I need to add ESRI image layer to map. So I now have my map's spatial reference as web mercarter. Now if I execute the GP Service nothing happens. No request being send to server, because nothing gets recorded in fiddler.

If I test the same request from ArcGIS Server rest end point it works fine.

Here is the request string in NAD27 and web mercarter. Only difference is wikid.

{"features":[{"geometry":{"type":"polygon","rings":[],"_ring":0,"spatialReference":{"wkid":4267,"latestWkid":4267}},"symbol":null,"attributes":{"id":1,"legalParams":"MT:TR|CNTY:051|STC:OH|PAR:|QQ:NW|LN:|PGUID:{3A15DFEB-9868-4588-9A08-737C033EDF2D}|MSC:1|STN:34|MER:01|TWSP:6|TDIR:N|RNG:5|RDIR:E|SEC:13|RW:|RUOM:AC"},"infoTemplate":null}],"spatialReference":{"wkid":4267,"latestWkid":4267},"geometryType":"esriGeometryPolygon","fields":[{"name":"Id","type":"esriFieldTypeOID"},{"name":"legalParams","type":"esriFieldTypeString"}]}

Web Mercator request

{"features":[{"geometry":{"type":"polygon","rings":[],"_ring":0,"spatialReference":{"wkid":102100,"latestWkid":3857}},"symbol":null,"attributes":{"id":1,"legalParams":"MT:TR|CNTY:051|STC:OH|PAR:|QQ:NW|LN:|PGUID:{3A15DFEB-9868-4588-9A08-737C033EDF2D}|MSC:1|STN:34|MER:01|TWSP:6|TDIR:N|RNG:5|RDIR:E|SEC:13|RW:|RUOM:AC"},"infoTemplate":null}],"spatialReference":{"wkid":102100,"latestWkid":3857},"geometryType":"esriGeometryPolygon","fields":[{"name":"Id","type":"esriFieldTypeOID"},{"name":"legalParams","type":"esriFieldTypeString"}]}


texts:
Is it possible to crop a .vrtmosaic to a vector-geometry? I have a .vrt-file which is a mosaic of various sentinel2-scenes. For me it seems to be more efficient to already crop and mask (set everything outside the polygon to NA) before converting to a .tiff-file.
I know that I can use something like gdalwarp -cutline... to crop the .tiff. But does this work with a .vrt aswell?


texts: How do I disable mouse scroll zoom in OpenLayers 3?

texts: I'm trying to override it 'cause it's so dang sensitive. Bonus points if you can also demonstrate a elegant way to override, but, well, I can probably figure that out.


texts:
I have a NASA BlueMarble tile layer that goes from zoom factors 0-8. I when I go to zoom level 9, I do see a momentary scaling to 9 (stretching pixels). But  then, as expected, I get the errorTileUrl.  Nothing wrong here.
But rather than have to created the tile cache (and massive data store) for zoom factors 9-14 or so, how would one go about having those zoom factors be pixel scaling and stretching?

This would be very useful for my offline maps.

Ideas?

(BTW, any suggestions for how to store offline tile caches? Say 2K to 4K of tiles in PNG format? Would I use an app manifest? (web application)?


texts: Use if statement to delete row data, ArcGIS 10.3.1, Python 2.7

texts:
I have a indoor positioning shp file. I need to do the data cleaning for the further analysis. and here is the attribute table:

Field(Sheet1_minor) is the name of beacon, and Fields(Beacon_41_B1,B2,B3,B4) which are near the beacon(Sheet_minor).
I am trying to clean the noise data by filtering the same record or adjacent record. For example, i start to work from ObjectID 28, I want to testify if ID 29 value of field(Sheet_minor) is equal to ID 28 adjacent record(Beacon47_B1,B2,B3,B4).

Yes->

testify ID 28 value of field(Sheet_minor) ==  ID 30 value of field(Sheet_minor), if this equation is true,delete ID29 and ID30. if this equation is false, skip it.

No->

skip this record.
Here is my code, but it doesn't seem to work.

# Licence:     <ArcGIS 10.3.1>
# Import python modules
import arcpy,sys
from arcpy import env

# set the workspace enviroment
env.workspace = r"D:\KKC\Indoor positioning\Ni_Feature\Ni_shp.gdb"
# To aviod an error, set the geoprocessing environment to allow existing data to be overwritten.
env.overwriteOutput = True
# Create a variable with the name of the subject feature class
fcName = "User1_test"
rows = arcpy.UpdateCursor(fcName)
fields = arcpy.ListFields(fcName)
# print out the fields in table
#    with arcpy.da.SearchCursor(fcName,["Sheet1__minor"]) as cursor:
for row in rows:
    for field in fields:
        if fields.getValue(field.Sheet1__minor) == fields.getValue(field.Beacon_47_B1) or fields.getValue(field.Sheet1__minor) == fields.getValue(field.Beacon_47_B2) or fields.getValue(field.Sheet1__minor) == fields.getValue(field.Beacon_47_B3) or fields.getValue(field.Sheet1__minor) == fields.getValue(field.Beacon_47_B4):
            if rows.getValue(field.Sheet1__minor) == rows.getValues(field.Sheet1__minor) + 2 :
                cursor.deleteRow()
            else:
                pass
                rows.updateRow(row)
        else:
            pass

-------------------------------------------
Traceback (most recent call last):
  File "<module1>", line 30, in <module>
AttributeError: 'list' object has no attribute 'getValue'
------------------------------------------

this is the message. I think i create the list for fields and rows.

Do you have any suggestion?


texts:
I want to use "Calculate field" based on the string in another field. I would like to loop if condition but something is wrong and I am little bit confused.

I am super new in Python. Here is my simple code:

import arcpy
stands = r'D:\mountaine\database.mdb\testPoly'
arcpy.MakeFeatureLayer_management(stands,'pormap')
with arcpy.da.SearchCursor("pormap", [ 'TEXTSTR', 'NUMBR']) as searchIng:
    for x in searchIng:
        if x[0] == 'BK':
            arcpy.CalculateField_management("pormap", "RESULTS", "[NUMBR]*1")
        elif x[0] == 'SM':
            arcpy.CalculateField_management("pormap", "RESULTS", "[NUMBR]*2")
        else:
            arcpy.CalculateField_management("pormap", "RESULTS", "[NUMBR]*3")

All I want is to multiply number from "NUMBR" field based on string in "TEXTSTR" field. Can someone tell me where I'm having a mistake ?

Here is my table:


texts: Cleaning up OpenStreetMap data, building polygons and trying to improve appearance

texts:
I'm currently going through my first project with OpenStreetMap data and I'm struggling to see how certain features can be made to look anywhere near as good as in OpenStreetMap.org.

I downloaded a pre-built Shp from Cloudmade for Ontario and many of the big lakes in the south are terrible - a polygon might cover a fraction of a lake's actual area and makes for poor visualisation. I think this is a result of OSM data only containing points and lines, and someone having to build polygons from these lines. Furthermore, I looked at the administrative boundaries lines layer and there were 100,000s of line features for all kinds of administrative entities but apparently no attributes that could help me build a single polygon for the province's full border.

Now, I know that the Cloudmade Shps have been through some post-processing to whittle the many feature types in the .osm XML down to 7 or 8 layers, so perhaps the way they collate administrative boundaries is not perfect for my area. Also, as the Ontario Shps don't cover the US states that also border the big lakes perhaps this explains why the lake polygons were not well formed, but if I have to download the .osm XML for such a large area and process it all myself this causes other issues - it will take a lot of time and many iterations and I don't know if the end result will be worth the effort (e.g. what if the lines that make up a lake's border don't properly meet up, and snapping introduces unwanted artifacts - or maybe none of the lakes are actually 'closed' because they all have rivers connected to them, so building a polygon is just impossible?)

I'm guessing that to make the map in OpenStreetMap.org someone must have processed all the data to build polygons where necessary. Presumably all country / state / province borders were created as polygons and shaded, and either the same was done for all lakes / oceans or the map simply has a blue background and everything else is layered on top? If this polygon data was created for the OpenStreetMap.org map, why can't I access it anywhere? It seems like there should be an additional resource between the raw XML and the state / province-specific Shps that have been processed according a particular rule-set

If anyone has come across similar issues and has any kind of advice I'd really appreciate it!


texts: hide and show a feature from group feature with checkbox in leaflet

texts:
To begin I will explain the code, I have a checkbox generated after loading the web page by the next code

    var CPID = new L.GeoJSON(false,
{
          onEachFeature: function (feature, my_Layer) {
           var element = document.getElementById("DDL");
           //alert(feature.properties.id_cp)
           element.innerHTML += '
&nbsp;&nbsp;&nbsp;&nbsp;<input type="checkbox" id="' + feature.properties.id_cp + '" value="' + feature.properties.id_cp + '" name="Circuit" onchange="CP(this.id)">';
           element.innerHTML += '&nbsp;&nbsp;' + feature.properties.id_cp + '';

}
});
var urlCPID = 'http://localhost:8080/geoserver/Gestion-proprete/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=Gestion-proprete:CP_ID&maxFeatures=500&outputFormat=text/javascript&format_options=callback:getJsoncp&srsName=EPSG:4326'
$.ajax({
    jsonp: false,
    url: urlCPID,
    type: 'GET',
    dataType: 'jsonp',
    jsonpCallback: 'getJsoncp',
    success: function (data) {
        CPID.addData(data);
    }
});

This code will fill the list of the checkbox from layer wfs share in geoserver after this

when I  checked one of this box he called this function

function CP(code) {
    var group = new L.featureGroup();
    var CPG = new L.GeoJSON();
    var CPG = new L.GeoJSON(false,
            {
                onEachFeature: function (feature, my_Layer) {
                    my_Layer.bindPopup("ID: " + feature.properties.id_cp + "
ID_arrondissement: " + feature.properties.id_a);
                }
            });

    var urlCP = 'http://localhost:8080/geoserver/Gestion-proprete/ows?service=WFS&version=1.0.0&viewparams=id_cp:' + code + '&request=GetFeature&typeName=Gestion-proprete:circuit_planifier&maxFeatures=1960&outputFormat=text/javascript&format_options=callback:getJson' + code + '&srsName=EPSG:4326';
 $.ajax({
        jsonp: false,
        url: urlCP,
        type: 'GET',
        dataType: 'jsonp',
        jsonpCallback: 'getJson' + code,
        success: function (data) {
            CPG.addData(data);
        }
 });

 if ($("#" + code).prop('checked')) {
     CPG.addTo(group);
     group.addTo(map);
 }
 else {
     map.removeLayer(group);
 }
}

In this code I have a query that take a parameter, id of feature then it make me to  display features by id. For the display it work very well but to display the feature checked its not working.
The problem is i want to display features by checking


texts: Join Attributes by Location incorrect summary attributes in joined table

texts: I am using QGIS 2.18.3 to join points to polygons with a Precision setting of 20.0 to select any polygon within radius of my point. I am taking attribute of intersecting features with the point layer as the target. The resulting joined layer attribute table confuses the field names. By checking results with a small test data set it is clear that the field names in the joined layer are displaced. I have 5 numeric fields in the polygon layer, two of which are Integer, the rest Real. The joined layer takes the summary values from the 3rd numeric field, but labels it as the 1st numeric field etc. I checked all values and apart from displacement all are correct. I am a novice user but is this a bug?


texts: Choosing aesthetic location over physical location?

texts:
I am in the process of putting together an infrastructure map for my city.
We have received some engineering maps (and CAD files) that are assisting me in my task.
I have also gone out and gathered GPS coordinate locations for some features (manhole covers, fire hydrants, etc).

Some features line up within a reasonable margin of error, but other features seem to be in locations fairly far off from where I found them in the field.
The features that are not in their "correct" locations seem to have been placed where they should be, in reference to other surrounding features.

An example of this would be a manhole cover that does not line up squarely with the intersection of two streets in the real world. In mapping software, things must connect. Therefore the manhole cover is aligned with the intersection of the streets for aesthetics and consistency.

For features like manhole covers and fire hydrants, I personally don't believe that we need pinpoint accuracy to determine where a feature is. For recognizable features such as those, who aren't usually clustered together, an accuracy of "it should be around here somewhere" seems accurate enough to me.

Is it a common practice to choose aesthetically consistent locations over physical locations?


texts: MSSQL performance with Geography (point within polygon)

texts:
I have a query, below, that works - however its very slow, looping through roughly 24,000 points per second. This is a rather high end production DB server with MSSQL and SDE, and in other instances a 10,000,000 point shape can be loaded in to view in ArcMap in ~5 seconds.

The question is: Can the SQL query be structured to perform any better?
My polygon is being given as points, and being converted to geography.
My data has a Shape column that according to MSSQL is stored as geography.
In order to get the query to work I'm asking every field to be converted via STAsText - and I think that's a performance hangup. For the life of me I cannot seem to get the query to work any other way. If I remove this, the query fails to run.

Ultimately a similar query will be run on a table with 20,000,00 points, so the performance as it is currently is not really acceptable (That would be 833 seconds).

The code as it is now:

DECLARE @geog1 geography;
SET @geog1 = (geography::STPolyFromText('POLYGON((153.013377 -27.385831, 153.018303 -27.399108, 153.030438 -27.395896, 153.028725 -27.386901, 153.013377  -27.385831))', 4326));

SELECT * FROM [GISDB].[GISUSER].[POINTS]
WHERE @geog1.STIntersects([POINTS].[Shape].STAsText()) = 1;

.

The data in the shape field is stored like such (Well Known Binary, I'm guessing)

SELECT [Shape] from [GISDB].[GISUSER].[POINTS]
WHERE [OBJECTID] = '100000';

    Returns: 0xBB100000010C300A689B1B4B3BC0400E645B061F6340


texts: How to apply "Band" settings using gdal Python bindings?

texts:
I am trying to replicate the gdal_translate command, but using pure Python:

  gdal_translate infile.tif
  outfile.tif -co COMPRESS=JPEG -co
  PHOTOMETRIC=YCBCR  -co TFW=YES -b 1 -b 2 -b 3

I have already looked at How to call gdal_translate from Python code? but am stuck with how to apply band information using the GDAL Python bindings.

So far I can see how to do all of this with the Python code, except apply the band (-b) settings.

#Import gdal
from osgeo import gdal

src_filename = r"infile.tif"
dst_filename = r"outfile.tif"

#Open existing dataset
src_ds = gdal.Open( src_filename )

#Open output format driver, see gdal_translate --formats for list
format = "GTiff"
driver = gdal.GetDriverByName( format )

#Output to new format
dst_ds = driver.CreateCopy( dst_filename, src_ds, 0, ['COMPRESS=JPEG','PHOTOMETRIC=YCBCR','TFW=YES'])

#Properly close the datasets to flush to disk
dst_ds = None
src_ds = None

After reading the GDAL API Tutorial, it appears that I can only do this with the Create() method, not the CreateCopy() method.  If I have to go the route of using the Create() method, what is the right approach to doing this?  It looks like using a VRT is likely the right way, but I cannot wrap my head around how to write the code that would "make a copy" of the existing file while using either the Create() or VRT methods.


texts: GeoServer 2.7.6 not displaying TIFF correctly

texts:
I would like to upload TIFF and display it as WMS layer. But after creating layer process, WMS layer is malformed. Please be lenient, I'm fresh in GeoServer. Maybe I did something wrong that I dont understand. Here is TIFF that I want to display as WMS:

https://www.dropbox.com/s/6dpoyvk372utjfd/Keihaskoski_dem_EDGE.tif?dl=0

WMS layer settings:

second:

third:

fourth:

malformed WMS:


texts: Layer vs. Color Transparency with ArcGIS Javascript API

texts:
When creating services with ArcGIS Server, I have paid attention to the warning that Layer Transparency will make response time slower and that I should change to color transparency instead so I have been changing it when I "share" the service. Maybe I'm not doing this correctly because I'm not getting the result that I want or maybe I need to ignore the warning.

Currently all my layers are in the same service and I can set opacity to some value 0.5 (and/or can set it individually for different layers e.g. just the layer on top) but the layers on top still cover layers on the bottom. I'd like to see the parcel boundaries underneath the other polygon layers.

mylayer[2].setOpacity(0.1);

Even if I set layer 2 to something almost invisible (like 0.1 as shown in code above; and as rendered below) it covers the parcels underneath (see parcel 0353 in the screenshot below). How should I be doing this instead? Do I need to ignore the warning?


texts: QGIS Exporting geometry of repeated labels of line layer?

texts: I would like to export geometry (in any format, preferably as a separate point shape, but will work around if necessary) OF LABELS of a line shape layer, where I use repeated labels ("Parallel" labels, repeated at a given "Repeat" distance). Is there a good way to do that, can anybody help?


texts: Find intersection point in a road network of line shapefile

texts:
How can I find intersection point in a road network in one line shapefile?

I am using QGIS 2.18.

How can I get the intersection point shown in blue colour


texts: Editing GML data in QGIS and also in ArcGIS?

texts:
I have a GML file with graphical and non-graphical layers, it opens in QGIS 2.11.14 but it doesn't allow anything but read permission. I can't edit new attributes or rewrite it... it seems like the GML is locked.

Why does it happen?

it also happens in ArcGIS...


texts: How can a HTML table be stored within a database field as a BLOB, text or other value?

texts:
My Spatialite database has several thousand spatial polygon records In QGIS, each with a couple dozen field values.  I also have many thousands more of additional tabular data relevant to those polygons, all stored within Spatialite, that the polygons need, a one to many relationship (a table as a value)

I'd like to create a QGIS ACTION of a single Field value that pops up a whole table stored within a single value in the spatial record.   My initial strategy is to use a BLOB or Text to store HTML coded tables I create, so that the QGIS ACTION would open it up in a web browser.  Opening up an equivalent table in Excel or some other common application would be satisfactory.

My worries are that I can't store html as a BLOB since its not binary.  (but maybe that's not relevant).  Also, I figure I can't store the html table as text since it would contain more characters of data than spatialite or shapefile can handle.

I can create an html coded page for each field value without a problem, but...

How would I store the html page as a single field value in Spatialite or Shapefile?
Or, what would be a good alternative to doing that, such as using another format besides html or BLOB data.
Would it be possible to have a SQL query stored as a field value, and assign an Action that pops open the Results?

I'm understanding I could use PostGIS to create a database relationship if all else fails.   It may be a day before I get back with any feedback.  Thanks in advance for everyone's suggestions.


texts: Create a field indicating if the polyline intersects the shapefile boundaries

texts:
I have a shapefile delineating 563 districts, I want to create a field that indicates whether the road network (red line) intersects the district or not.

How do I do this?


texts: How do I keep layer scale visibility when exporting a leaflet map with qgis2leaf?

texts: I'm trying to export a leaflet map using qgid2leaf. Some of the map layers have a scaled visibility range. However on the leaflet map these layers show at all scales? Can this be done?


texts: Calculate Width and Height of a Google Static Maps Image

texts:
Do you know how to calculate the width and height (in degrees minutes) of a Google Static Maps image?

I've found it extremely difficult. The only information you have is:
- The centre of the image (as a Longitude Latitude coordinate)
- The zoom level of the satellite (1 >= zoom <= 25)
- The width and height of the image in pixels

See this image for a visual explanation of the problem:

And from this limited information its very difficult to calculate the image dimensions in degrees minutes, mainly because the zoom level heights (in feet or metres) is not exposed by Google AFAIK. AND because you can't move up on the y axis at a constant value because you have to take into account that the earth is a sphere and not a flat surface so moving up at a constant value won't work (which is what all the answers on StackOverflow provide). It appears meridional points is a solution but the math involved is completely beyond me.

Any advice would be greatly appreciated.


texts: Changing CRS using GeoPandas in Python when print(file.crs) gives 'None'

texts:
I'm trying to change the .crs from a cylindrical projection (WGS84 (lat/lon)) to a Mercator-projection.

Some information can be found here. However it doesn't seem to work for me for this shapefile of Belgium. (the example on the GeoPandas website for the world worked well so all libraries are installed correctly)

Someone an idea what the problem might be? -> My .crs stays cylindrical and does not change to Mercator-projection for this shapefile of Belgium.

Dataset: 'BELGIUM__Municipalities'

Example-code:

import geopandas
import fiona
import matplotlib.pyplot as plt
import pandas as pd

def records(filename, list):
    list = sorted(list)
    with fiona.open(filename) as source:
        for i, feature in enumerate(sourceô:max(list)+1):
            if i in list:
                yield feature

#a = list(range(588)) #Belgium
a = list(range(70)) + list(range(89,154)) + list(range(181,310)) + list(range(463,507)) #region of Flanders in Belgium

municipalities = geopandas.GeoDataFrame.from_features(records("BELGIUM__Municipalities.shp",a))
print(municipalities.crs) #has None as output! -> so I have to set a crs myself
municipalities.crs = "epsg:4326"  #WGS84(lat/lon)-projection
municipalities.plot(facecolor = 'lightgrey', linewidth = 0.05, edgecolor = 'black', alpha = 0.25)

municipalities = municipalities.to_crs("epsg:3395") #Mercator-projection
municipalities.plot(facecolor = 'lightgrey', linewidth = 0.05, edgecolor = 'black', alpha = 0.25)

plt.show()


texts: How to define correctly the Grass region in Qgis

texts:
I know that my question could be pretty basic for some of you. I am trying to learn GRASS. I have some raster layers and a vector file in  EPSG 32614. Using the grass plugin I  defined a new map set, defined a region using the limits of my raster layers.

I succesfully imported both the vectors and the raster in to grass. When I try to run r. cost I have the following error:

Raster map is in different zone (14) than current region (0)

￼Unable to open raster map

Nevertheless, the metadada of my raster and vector layers are in  32614 (UTM zone 14N). The grass region settings tool tell me this :

I looked in to this question

Raster map in different zone () than region ()

But It does not tell how ti fix the problem. I would appreciate any hint on the matter. I have the same problem working with 2.1 64 bits standalone and also in Lisboa 32 bits standalone working under win 7 64 bits.

Thanks in advance

Gerardo


texts: Creating CRS with coordinate origin at north-west corner of Mercator projection?

texts:
I want to define a custom CRS for a very special, technical use-case. For that I use a WGS84-World-Mercator as base definition, here are the proj4 parameters taken von QGIS:

+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

For my custom CRS, I wish to have positive coordinates for any point in the world and all coordinates should be within the range of 0 and 1. The coordinate origin should be in the north-west corner of the Mercator.

To scale all coordinates within the range [0, 1], I use the parameter +to_meter=40075016.68557849 while leaving the scale factor k_0=1 and removing the metric units=m.

+proj=merc +lon_0=0 +k_0=1 +to_meter=40075016.68557849 +x_0=0 +y_0=0 +datum=WGS84 +no_defs

So far, so good, now the tricky part: How to create a coordinate origin at the north-west corner (North Pole)? Taking pen and paper, I found the solution is setting +lat_0=90 and lon_0=-180, like that:

+proj=merc +lon_0=-180 +lat_0=90 +k_0=1 +to_meter=40075016.68557849 +x_0=0 +y_0=0 +datum=WGS84 +no_defs

However, that creates weird results:

The +lat_0=90 basicly does nothing to whatever value it gets set. Why?
The +lon_0=-180 does behave like setting it to plus +180°, but I explicitly want the negative value from a mathematical view.

I know that plus and minus 180° are the same meridian on the map, but technically the difference with the minus is important for me.

How to create a CRS with coordinate origin at the north-west corner of the Mercator projection and only positive coordinates throught out the whole map? Is inverting the y-axis the issue here? Is chosing a World Mercator as base the wrong approach?


texts: Generating bounding box based on tile numbers using slippy map tiles method

texts:
Description:
I am writing a map cache server for my company. I have to get input data based on tms standard and send a request to wms server to render and cache map files.
I am using slippy map tiles formula to convert x,y,z to latlon. (North West)
Also for finding the other corner (south East), I am converting x+1 and y+1 to degrees. We are using srs 4326

My problem:
It works basically, but the map files are a little scaled up! in lower zoom levels, the problem is very obvious.


texts:
I have been trying to create a torque animation where the points on the map fade out at different times.  I have a table with 'start_date', 'end_date', and a 'day_count' which is a count of the number of days between the start and end date.

Is there a method of setting up the Torque visualization so that when the time slider gets to the end dates for each point the points will disappear. I have not seen any similar examples on cartodb.com, or github.


texts: ArcPy: Excel to Table Conversion and Join Management

texts:
I am trying to join an Excel table to a shapefile in ArcPy. My first step is to convert the Excel sheet to a table, but I am receiving an error:

#Convert xlsx to table:
import os
import xlrd
import arcpy

in_excel= r'C:\\Temp\\Planning\\Planning_Out\\JoinTest.xlsx'
out_gdb= r'C:\\Temp\\Planning\\Planning_Out\\JoinTest.gdb'

def importallsheets(in_excel, out_gdb):
    workbook = xlrd.open_workbook(in_excel)
    sheets = [sheet.name for sheet in workbook.sheets()]

out_table=  os.path.join(out_gdb,arcpy.ValidateTableName("{0}_{1}".format(os.path.basename(in_excel), sheet),out_gdb))

print('Converting {} to {}'.format(sheet, out_table))

        # Perform the conversion
arcpy.ExcelToTable_conversion(in_excel, out_table, sheet)

Error:

Traceback (most recent call last):
  File "C:\Temp\Planning\LAMP.py", line 52, in <module>
    out_table=  os.path.join(out_gdb,arcpy.ValidateTableName("{0}_{1}".format(os.path.basename(in_excel), sheet),out_gdb))
NameError: name 'sheet' is not defined

My code was taken from:
http://pro.arcgis.com/en/pro-app/arcpy/get-started/working-with-geodatabases.htm

My thought process is that once the Excel file is converted to a table, I can join it to my shapefile using the following code:

# Make the shapefile a layer
arcpy.MakeFeatureLayer_management("merged_lines.shp", "tempLayer")

# Set the local parameters
inFeatures = "tempLayer"
joinField = "Name"
joinTable = r'C:\\Temp\\Planning\\Planning_Out\\JoinTest.gdb '
fieldList = ["Name", "Status"]

# Join feature classes by joinField
arcpy.JoinField_management (inFeatures, joinField, joinTable, joinField, fieldList)

I am using version 10.3.1


texts: gdaldem - raster to RGB raster with color-relief

texts:
I have 6-layer raster with float value.
I have six different QGIS style with color value.
Here is image value from QGIS:

Here is a QGIS style (for 6 layer):
<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
<qgis version="2.18.20" minimumScale="inf" maximumScale="inf" hasScaleBasedVisibilityFlag="0">
  <pipe>
    <rasterrenderer opacity="1" alphaBand="0" classificationMax="inf" classificationMinMaxOrigin="Unknown" band="6" classificationMin="-999" type="singlebandpseudocolor">
      <rasterTransparency/>
      <rastershader>
        <colorrampshader colorRampType="DISCRETE" clip="0">
          <item alpha="255" value="-999" label="-999" color="#a9a9a9"/>
          <item alpha="255" value="5" label="0-5" color="#0000ff"/>
          <item alpha="255" value="10" label="5-10" color="#00f400"/>
          <item alpha="255" value="20" label="10-20" color="#ffff00"/>
          <item alpha="255" value="30" label="20-30" color="#e0a000"/>
          <item alpha="255" value="inf" label=">30" color="#ff0000"/>
        </colorrampshader>
      </rastershader>
    </rasterrenderer>
    <brightnesscontrast brightness="-12" contrast="0"/>
    <huesaturation colorizeGreen="128" colorizeOn="0" colorizeRed="255" colorizeBlue="128" grayscaleMode="0" saturation="0" colorizeStrength="100"/>
    <rasterresampler maxOversampling="2"/>
  </pipe>
  <blendMode>0</blendMode>
</qgis>

I convert it to color-relief file
-999 169 169 169 255
0 0 0 255 255
5 0 244 0 255
10 255 255 0 255
20 224 160 0 255
30 255 0 0 255

I use gdaldem:
gdaldem color-relief -b 6 input.tif color-relief.txt output.tif

After I run this is problem with interpolation:
Correct image in QGIS:

Incorrect image from gdaldem:

How to set legend with interval value and disable interpolation of pixels?
I use -nearest_color_entry flag. It works but compartment is wrong.
gdaldem color-relief -nearest_color_entry -b 6 input.tif color-relief.txt output.tif
-999: 169 169 169 255
0-5: 0 0 255 255
5-10: 0 244 0 255
10-20: 255 255 0 255
20-30: 224 160 0 255
30-9999: 255 0 0 255


texts: Importing shapefile (.shp) in Adobe Illustrator

texts:
I am working in ArcGIS, I need to import the shapefiles with all layers, to adobe illustrator.

Can anyone suggest the steps doing this?


texts: leaflet.heat heat map highlighting overlapping areas

texts:
Lately I've been working on heat maps using the leaflet.heat plugin.

I have a set of points in the WGS84-system with the following structure:

latitude | longitude | weighting | type

The weighting-attribute is currently set to 1 so it will not influence the occurence of the heat map.

The type-attribute defines, by an number, which type of animal the dataset is about, let's say 1 stands for zebras and 2 for lions.

There is no problem in displaying these two heat maps seperately (as in the first two pictures), however when I'm trying to display both maps simultaneously it doesn't fit my expectations (third picture).

What I'm interested in is the area where those two heat maps overlay, namely in the middle of northern and southern lakeside, this is the area that should be highlighted the most.

Is there any way to do this with this plugin?


texts: Add Polygons to SpatialPolygons via loop / iteration using R?

texts:
Just out of interest: I'm creating a SpatialPolygon (R with package "sp") with a number of Polygons within. Is there a way to loop this task? I'm doing it manually as follows. The Polygons consist of 5 points each (i.e. a square) that are stored together in a larger matrix ("bbox_matrix_sp").

SpatialPolygons(list(
    Polygons(list(Polygon(bbox_matrix_sp[1:5,])), "reh1"),
    Polygons(list(Polygon(bbox_matrix_sp[6:10,])), "reh2"),
    Polygons(list(Polygon(bbox_matrix_sp[11:15,])), "reh3"),
    Polygons(list(Polygon(bbox_matrix_sp[16:20,])), "reh4"),
    Polygons(list(Polygon(bbox_matrix_sp[21:25,])), "reh5"),
    Polygons(list(Polygon(bbox_matrix_sp[26:30,])), "reh6"),
    Polygons(list(Polygon(bbox_matrix_sp[31:35,])), "reh7"),
    Polygons(list(Polygon(bbox_matrix_sp[36:40,])), "reh8"),
    Polygons(list(Polygon(bbox_matrix_sp[41:45,])), "reh9"),
    Polygons(list(Polygon(bbox_matrix_sp[46:50,])), "reh10"),
    Polygons(list(Polygon(bbox_matrix_sp[51:55,])), "reh11"),
    Polygons(list(Polygon(bbox_matrix_sp[56:60,])), "reh12"),
    Polygons(list(Polygon(bbox_matrix_sp[61:65,])), "reh13")
    )) -> sp_re_alle

This seems a little bit tedious, especially if the number of Polygons within the SpatialPolygon increases.

Is there a faster, more elegant way?


texts: Creating additional attributes for regional intersections?

texts:
A certain country is (simultaneously and independently) divided into 400 regions and 100 judicial districts. Each division is represented by a separate multipolygon layers with units as features. Intuitively one judicial district may contain many regions. However, one region can belong to many judicial districts (eg. Northern part of the region R_1 can belong to the district D_1 but its southern part can belong to the district D_2) - the aim of this is to reduce the risk of collusion between local judges and authorities.

I would like to create a database describing the relations between regions and districts. In the “districts” layer I’d like to create 400 columns/attributes representing all regions. The values of the attributes should be binary i.e. “1” if a given region overlaps/intersects with a given district and “0” if it does not.

For instance, for the district D_37 I would like to have “O” value of the attribute R_123 because the “region 123” does not overlap with the district 37 and “1” for the attribute R_44 because the “region 44” does actually overlap with the district 37.

How I should do this to avoid unnecessary difficulties (I use QGIS)?


texts: How to remove white stripes in TIN generated in Arcgis?

texts: I've generated TIN from DEM data using toolbox (3Danalyst/conversion/from raster/raster to tin) in arc-desktop 10.2. I'm getting these stripes, I tried altering toolbox parameter values but couldn't solve it, could anybody help?


texts: count if a row value is within a dynamic range?

texts:
I am after doing a COUNTIFS in python. I need to check each row ± 250 with other rows and if this row is between that range then count 1, then next row, .. and do the same thing. and if there are any, put the total number, else put 1. Here's the code in excel: COUNTIFS(D:D,D2,E:E,">="&E2-$G$2,E:E,"<="&E2+$G$2)

Here's the code, but it's not doing countifs.

import arcpy
fc = "feature_class"
counter = 0

with arcpy.da.UpdateCursor(fc,("Mean_m","Count")) as cursor:
    for row in cursor:
        if row[0] >= row[0]-250 and row[0] <= row[0]+250:
            counter += 1
            row[1] = counter
        else:
            counter = 1
            row[1] = counter
        cursor.updateRow(row)


texts: QGIS GEarthView plugin ReactorAlreadyInstalledError

texts:
I've been battling to get this really handy plugin working for some time. When I try installing it from plugin manager, I get the following error:

  "This plugin is broken
  'module' object has no attribute 'ReactorAlreadyInstalledError'"

I've just installed QGIS 2.8.1 64Bit on Windows 7 machine. Zope and Twisted libraries both installed.

Any ideas?


texts: The geotools library's GeoJSON plugin requires an ID to generate a Feature type. The specification of GeoJSON says it is a SHOULD and not a must. I'm trying to get rows from a database like SQL Server or Postgres and convert it to GeoJSON. Should I make the ID to the primary key of the table or I should let generate the ID by the library?


texts: I have a linear network of lines all digitized in the direction of flow and am trying to assign route ID's to segment groups. Is there an automated process to do this using QGIS?


texts: Where can I download a xyz bathymetry file of the Northwest Atlantic?

texts:
I'm looking to download the bathymetry of the Northwest Atlantic ocean as a txt file or any other file format that I can easily upload and use in R.

I once used the Global topography - Scripps Institution of Oceanography but it looks like I'm running into an error message this time...?

I have been looking at other websites, but I'm having a hard time to locate a new site where I could download a xyz file of my area of interest.

I would appreciate any recommendations.
Thank you very much in advance,


texts: How to get all the Rules and Filters from GeoServer for specific Workspace

texts:
I am using GeoServer 2.4.2 for storing my vector data. And I created a style for a layer
now I want to fetch all the Rules and Filters from style which I created on my work-space.
I am currently using this url for getting output :

http://mydomain.com/geoserver/rest/sldservice/{layerName}

But I am getting a default screen as shown :

Can anybody please let me know is there any way to get all Rules and Filters from GeoServer for specific Workspace ?


texts: Countries' Latitude and Longitude range?

texts:
Do countries have a Latitude and Longitude range?.

Is it possible to get the latitude and longitude range or boundary for each country.

UK lat and long is 51.5000° N, 0.1167° W.

What is the Maximum and minimum ,latitude and longitude i can have to be still in the UK?

For an example, i have 2 sets(a,b) of lat and long, I start at set a which i know is within the country but i don't know if Set b my destination is still within the country.


texts: Is it possible to set different scale dependent visibility for one layer classifications?

texts: I'm working on a sewer and water network, so I need to set different scale visibility for primary, secondary and tertiary network.


texts: Creating polygon features composed of scalable points

texts:
I want to represent a forest using symbols I create (stylised trees).
I considered two options:
- A point feature.
- A polygon feature.
This is the one I'd use, as a forest is better represented with an area rather than with a cluster of points. Among the ESRI style symbols there are some which are close to what I'm looking for (e.g. "Scattered Trees 1" symbology).
I want to create a similar symbology using a picture I created, like a small tree (and have it multiplied homogeneously across the polygon area).
I know that it is impossible to scale such symbols (the small trees composing the polygon filling) when zooming in and out. But Point features can be scaled. Is there a way to distribute point features in order to cover a specific area (as if I were drawing a polygon)?


texts: Casting LineString to Polygon results in ClassCastException when y coordinates are the same

texts:
final WKTReader wktReader = new WKTReader();
String geomString = "MULTILINESTRING ((-98.753985 34.63093, -98.7555568 34.63093))";
MultiLineString geom = (MultiLineString)wktReader.read(geomString);
LineString newGeom = new GeometryFactory().createLineString(geom.getCoordinates());
Polygon envelope = (Polygon) newGeom.getEnvelope();

Such code results with ClassCastException(). However, when I'm changing y value of any of coordinates, e.g. first one to 34.63094, I don't get any error.
I am using org.locationtech.jts-1.17.
What would I like to achieve is to cast it to Polygon without this Exception.
Or, the exception is actually valid and I shouldn't be able to cast any two-coordinate LineString to Polygon?


texts:
I am completely new to QGIS (and mapping in general), and I can't get my layers to line up.

  QGIS Version: 3.2.1-Bonn Standalone-installer Windows 10 64-bit

All of my layers are in the USA Contiguous Albers Equal Area projection. Currently, I have three layers:

Site Coordinates
World Imagery layer from ESRI
Cropland Data Layer (CDL) from USDA-NASS CropScape

My coordinates line up perfectly with the World Imagery layer. However, the CDL is floating somewhere out in the universe and won't line up.

  CDL Extent:
  48885.0000000000000000,2818725.0000000000000000 :
  117585.0000000000000000,2930685.0000000000000000

I downloaded the CDL as a zipped file containing a .tif and a .vat.dbf.tif file.
I unzipped the file and uploaded it to QGIS as a .tif file.

Is it possible that the georeferencing for the .tif file is encoded in the .vat.dbf.tif file? How do I import the Cropland Data Layer to QGIS?


texts: Set up tool to only accept the name of a specific feature class?

texts:
I'd like the red x to show in the corner of the tool if they try to load the wrong feature class. That way we can prevent possible errors with the output.  I've tried using the error messages arcpy functions but I only get the error message function after I select ok for the tool to run.

def getParameterInfo(self):
# First parameter
in_features = arcpy.Parameter(displayName="Data",name="Data",datatype="DEFeatureClass",parameterType="Required",direction="Input")
messages.addMessage("Must Use Data Feature Class!!!!!")
clearMessage()
return [param0]

This will set the parameters so that any other feature class that isn't named "Data" can not be used.  However, it gives the following error message.  I would prefer to be able to customize the red x message saying that you can only input the "Data" feature class.

setErrorMessage

The Set Error Message also gives the invalid drop item warning.  I'd rather have it as the red x in the corner and the pop up box.


texts: How do I rotate the OSM map before printing?

texts:
What is the simplest way to rotate and render Openstreetmap data, and still get a nice render? For example Toronto, would fit best rotated slightly on a page in landscape mode, so the northern boundary of the city are parallel to the top edge of the printed page.

Example: Toronto rotated here with Gimp, works fine but the labels are not really supposed to be rotated.

The most common thing for Openstreetmap is to import the data into  PostgreSQL+postgis and render with Mapnik. But really any solution is fine, e.g. doing it in ESRI tools.

PS. This is a repost of a question asked on help.openstreetmap.org, I'll repost the answers there if I'm allowed.


texts: Why does the location of the mean of a set of points not coincide with the location of a point I created in QGIS?

texts:
So I've made a number points for a set of locations in the Netherlands. Here they are:

The location of the purple point was computed by means of the “mean coordinate(s)” tool in QGIS. It represents the mean of the coordinates of five student cities in the West of the Netherlands: Amsterdam, Leiden, Delft, Rotterdam, and Utrecht.
I wanted to create a point at the exact same location of this mean, and add it to the existing table of (locations of) student cities I already had. So I looked up the locations of the mean:

And I created a point for which the longitude matches the MEAN_X and the latitude matches the MEAN_Y. Here it is as point number thirteen in the table:

On the map, it corresponds to the yellow (highlighted) location. However, as one can see, the locations of the two points don't match. The points do not overlap one another.
My questions are: why is this the case, and how can I adjust the situation in such a way that the points do overlap? It seems that the coordinate systems of the two points match.


texts: QGIS Map to Autocad (not SVG,PDF OR Image)

texts: Is there any way to export the map extent set with Altas points in QGIS to convert to DWG instead of whole project or any shapefile?


texts: How to find length attribute form my all lines in my layer

texts:
I am a begginer for QGIS
I have a column for the length of the lines, but don't know how to find it. Can anyone help me?


texts: Using Iterators in ArcGIS Pro ModelBuilder?

texts:
I am currently struggling with two issues in ModelBuilder both dealing with iterators. First I have situation where I have two model (Model1 and Model2). Both of these models produce a series of raster outputs (Model1_A, Model1_B, Model1_C...) and (Model2_A, Model2_B, Model2_C...). I need to sum the outputs from each model together so that Model1_A get added to Model2_A, Model1_B to Model2_B ect. as paired inputs. Is there a way to do this using iterators in ArcGIS Pro (Esri removed the run variable as list option from ArcGIS Pro). Additionally how would I create a workflow where every possible combination of outputs from Model1 and Model2 are summed?

!enter image description here


texts: Conditional Labeling based on attributes in ArcGIS 10

texts:
I am trying to label a street centerline layer with a condition based on street width.  I am using the Python parser and I can't figure out this simple logic.  Neither of these questions address how to accomplish this within the "Label Expression" window of ArcGis using Python:

Customizing label features using arcpy

Displaying combined data values in ArcMap

My attempt at it looks like this:

def FindLabel ( [NAME] ):
  if [WIDTH] > 30
  return [NAME]

GIS returns me an error when I try this.

Any help is much appreciated.
Thanks in advance,
mike


texts: I downloaded some census tract data which are stored in .gdb. I am able to open them in ArcGIS as shpfiles. I selected some tracts and tried to export them into another folder. But I kept getting this error: "Error exporting data. The maximum record length has been exceeded". Is this because the attribute table has too many columns? What is the solution?


texts: Unable to execute ArcPy module outside of ArcGIS environment

texts:
I am using ArcGIS for Desktop 10.1 (Arc Info license) and ArcGIS for Server 10.1 and automating few things on server level and getting interesting errors.

I am overwriting map service through this script. The scripts works fine in ArcGIS Desktop environment. When I tried to run script outside of ArcGIS i.e cmd, it is getting an error.

I tried following things at my end:

My environmental variable is set properly
I have imported all ArcPy, OS modules & toolbox module.
After getting an error service is not overwriting.
Researching on my side.
Script code works fine in the ArcMap

Still getting an error & wondering that it is working greatly in ArcMap but outside of ArcGIS why it is not working?

Any help will be great !! thanks in advance : )

Please check complete code and screenshots


texts:
I could use some tutorials on remote sensing, but some books provided they are free, would be useful.

To be specific, tutorials on ERDAS Imagine, ENVI and other remote sensing software, general tutorials, for beginners but for intermediate and advanced levels as well.


texts: Is it possible to install PyQGIS in virtual environment? I want to use rendering capabilities in my own application which runs using virtualenv.


texts: I have a database in PostgreSQL / PostGIS with polygons so how can i merge some polygons into multipolygons of the database using it as QGIS Layer?


texts:
I just ran the OSGeo4W installer to update my previous installation and I noticed that it was updating my QGIS installation from 1.8 to 2.0.  I haven't been able to find any information on the web about QGIS 2.0 being released yet.

The installer ran without any issues, but now that it is complete I cannot find executables for QGIS.  The folder for QGIS exists within the apps folder under C:/OSGeo4W/, but there are no executables in there.

I've already completely uninstalled all OSGeo4W applications and reran the installer using the Express Install option, but I had the same result.

How can I use QGIS and has version 2.0 been released yet?


texts: Average distances between two point sets bound by polygon

texts:
I have a table of boundary polygons, a table of points and a table of other_points

I want to work out the average distance between points contained within and other_points for each polygon.

So far I have this query which works but is not very efficient, I limit distances to within 10,000 meters (data is in SRID 27700). I am using POSTGIS to do this, I ultimately need to get the data into Postgresql.

I am sure it can be done better as this will take many hours to compute. Any ideas how?

WITH POINTS AS (

    SELECT  point_id,
            polygon_id,
            point_geometry as wkb_geometry
   FROM boundary
   LEFT JOIN points

  ), DISTANCES AS (

    SELECT
        distinct ON (ogc_fid) ogc_fid,
        polygon_id,
        round(ST_Distance(B.wkb_geometry, D.wkb_geometry)::NUMERIC,1) as distance
    FROM
       POINTS B, other_points  D
    WHERE ST_DWithin(B.wkb_geometry, D.wkb_geometry,10000)
    ORDER BY ogc_fid,polygon_id,ST_Distance(B.wkb_geometry, D.wkb_geometry)

)
SELECT distinct polygon_id,
       count(polygon_id) over(partition by polygon_id),
       round(avg(distance) over(partition by polygon_id)::NUMERIC, 1),
       FROM DISTANCES
   WHERE DISTANCES IS NOT NULL


texts: I have a shapefile of buffers (n=10) and a shapefile of watersheds (n=10). Both shapefile contain overlapping features. I would like to clip each watershed feature with it's corresponding buffer, based on a common ID field (highlighted area below). I'm thinking ModelBuilder should be able to do this, but haven't been able to figure it out. I also have each watershed as an individual file if that helps.


texts: Assign different values to each zone, ArcGIS 10

texts:
Is it possible to assign different values to each zone. The value would be decimal.

100==> 4.97
200==> 5.90
300==> 6.01
..ect...
1000==> 9.64
or somehow change value with cost raster.
I have DEM raster, buffer raster, reclass raster, poligon layer.


texts: ArcGIS tool to batch re-order directions of polylines

texts:
I'm working with multiple datasets that each have thousands of polylines in their respective shapefiles. Is there an ArcGIS tool that allows me to batch re-order the directions of the lines? For example, order all lines N to S, then W to E. Will ultimately do this in ArcPy, so I don't want to use editor.

Example screenshot below. The lines were digitized in different directions (sometimes "top-down" and sometimes "bottom-up"). I just want to set a uniform rule.


texts: How to set up Geoserver app-schema feature chaining?

texts:
I am using Geoserver 2.12 with the application schema plugin to map a PostGIS database schema to the GeodesyML schema. I need to set up a WFS request that returns GeodesyML output such as:

<wfs:member>
  <geo:GeodesyML gml:id="geoml_monuments.fid-39d8a647_15cd161ff6e_-7e89">
     <geo:Site>
        <geo:installedDate>2001-09-12T06:00:00Z</geo:installedDate>
     </geo:Site>
     <geo:Monument>
        <geo:height />
        <geo:foundation />
     </geo:Monument>
     <geo:siteLog>
        <geo:siteIdentification>
           <geo:siteName>Cordell, OK</geo:siteName>
           <geo:fourCharacterID>SG19</geo:fourCharacterID>
        </geo:siteIdentification>
        <geo:siteLocation>
           <geo:city>Norman</geo:city>
           <geo:state>Oklahoma</geo:state>
           <geo:countryCodeISO>United States</geo:countryCodeISO>
           <geo:tectonicPlate>North America</geo:tectonicPlate>
           <geo:approximatePositionITRF>
              <geo:cartesianPosition>
                 <gml:Point gml:id="cart.geoml_monuments.fid-39d8a647_15cd161ff6e_-7e89" />
              </geo:cartesianPosition>
              <geo:geodeticPosition>
                 <gml:Point gml:id="latlon.geoml_monuments.fid-39d8a647_15cd161ff6e_-7e89" srsDimension="2">
                    <gml:pos>-98.977854677 35.355545782</gml:pos>
                 </gml:Point>
              </geo:geodeticPosition>
           </geo:approximatePositionITRF>
        </geo:siteLocation>
     </geo:siteLog>
  </geo:GeodesyML>

I have succeeded in doing this when the data is all stored in a single materialized view in my database. Unfortunately this is too limiting and I need to be able to pull things from multiple views. Specifically I want to be able to pull the geodeticPosition and cartesianPosition geometry objects from a PostGIS view called location_view and all the rest of the data from site_view. I believe the way to do this is by Feature Chaining wherein I need to set up a top level mapping file which pulls data from site_view and <Include>s a second mapping file that pulls from location_view. In each of these mapping files, what do I need to have in the <targetElement>, <targetAttribute>, <linkElement>, <linkField>, and <sourceExpression> fields? Alternatively, is there an easier way to pull data from two different tables/views into a single mapping file?

Here is an example provided from the Geoserver docs on Feature Chaining


texts: Postgis: find geometry within multiple polygons

texts:
I've got two tables with geometry in it.
Table 1:items holds items with a point geometry; while the second table property_zones holds polygons.

The Items:

CREATE TABLE public.items
(
  id integer NOT NULL DEFAULT nextval('item_id_seq'::regclass),
  ....
  geom geometry(Point,4326),
  ....
)

The zones

CREATE TABLE public.property_zones
(
  id integer NOT NULL DEFAULT nextval('property_zones_id_seq'::regclass),
  prop_id integer NOT NULL,
  ....
  geom geometry(Polygon,4326)
)

Both tables have GIST indexes on the geom column.

I now want to query these tables and find every item within the zones op a particular property (denoted with prop_id). I came up with this SQL, which works but has a flaw in it:

SELECT i.id
FROM items AS i
WHERE
   ST_Within(i.geom ,
               ( SELECT geom FROM property_zones WHERE prop_id = :propid )
             )

This works perfectly, but some properties have multiple zones, and off course the query within the ST_Within delivers more than 1 tuple. And of course I can return to a php solution where I split the query up in a couple of query's; but I really want a query based solution. (If possible)

How can I rewrite this code to find every item within the zones (1 or more) of each property?


texts: How to get the source name of a feature layer

texts:
I have feature layers that have a different name in the TOC then that they do in the GDB. I cannot change either of these and I need the name of the feature class from the GDB in order to make my model work. I have tried things such as parse path but I only get the name of the layer as it is in the TOC.

I have tried a few ArcPy/Python methods, but with my little knowledge of these things, I could not figure out how to get the GDB name as an output parameter.

I am using ArcMap version 10.3.1

Below I have added an example of my model and the code I last used to try and get the name.

import os
import arcpy

layer = arcpy.GetParameter(0)

print(os.path.basename(layer))


texts:
I'm having a problem with ArcGIS Server local edits. We have a feature service we use for multi-user editing and we use ArcMap 10.4 for creating a local copy and synchronizing when we are done.

Problem is, as we all know, sometimes ArcMap crashes. This would not be a problem if the user would save the MXD so they can sync back when they are done but sometimes this is not done and we are met with a rather difficult situation.

A simple and rather raw way to fix this is to find the local copy in FeatureServiceLocalEdits and replace the data on the server with the current data.

This can lead to rather frustrating tasks and potential data loss.

Is there a way to sync a local copy to ArcGIS Server after the group is lost?


texts: Batch import of DEM tiles (with shapefiles) into PostGIS creates non continuous contour lines

texts:
I want to import a large amount of ASTER DEM data into PostGIS (covering Europe). I use this script: https://wiki.openstreetmap.org/wiki/Contour_relief_maps_using_mapnik#Importing_ASTER_GDEMv2_data

The script creates the contour lines (gdal_contour) for every DEM tile in the form of a shapefile, then appends the data (shp2pgsql) to the contour table in postgis. The only thing I added is -s 4326:900913 to the shp2pgsql command.

Everything works, except the contour lines are not continuous at the tile borders (the borders of the DEM tiles):

The DEM tiles overlap by one pixel. I also tried to crop it, then there is no overlap like in the picture, but an offset between the contour lines.

How is it possible to import the DEM data and get continuous contour lines at the borders of the DEM tiles?

PS: I tried to merge all DEM tiles and generate one big shapefile. But gdal_contour aborts at the file size of 4GB.


texts: Convert vector data with irregular points (no grid) to raster with QGIS?

texts:
I am aware that how to convert a vector layer to a raster layer is a quite frequent problem and that there are many questions and answers dealing with it. However, none was helpful concerning my specific problem, so I hope it's okay to ask again.
I have a huge dataset consisting of fish catch data - different species and the location where they were caught. The original data is a csv file, but I'm working with it in QGIS and have created a vector file from the spreadsheet.
Now I want to use another program (Zonation) that only accepts raster files to calculate species richness in certain areas - is there a possibility to create a raster from a vector layer with irregular points? I have tried the Rasterize tool, but it only gives me a solid rectangle and all individual data points are lost. I'm working with CRS WGS 84 and I have tried to set the cell size to 0.4, but still it's only one solid square.


texts: Capabilities of the Analytical Hierarchy Process (AHP) in conjunction with weighted overlay tool

texts:
I am making a prospective map of fossil deposits using the "Weighted Overlay" ArcGIS tool. Input factor rasters are mostly topographic in nature (elevation, roughness, TPI, etc). Initially, I  had planned on assigning the weightings of each factor based on expert opinion. However, due to the large number of input rasters it is now apparent that a less subjective, statistical technique of determining these weights would be preferable. I want to know if the analytical hierarchy process (AHP) can be suitably applied to determine the weights of the input rasters in this case.

If not, can you suggest any other methods to reduce the subjectivity of the designation of these weightings?

I am using ArcMap 10.2


texts: Shifting HEC RAS cross sections to the side

texts:
I was wondering if there is a way to shift a cross section in HEC RAS to the side?
Since my data is not the best, a lot of xs are shifted somehow (see screenshot). The only way fixing it I found so far is opening and editing (shifting) the entire xs in RASMapper, then going back to the xs editor, checking whether the shifting is okay and depending on whether it does not look fine redoing it several times.
Would be okay when only having a bunch of xs, but the model includes hundreds of xs...
Does anyone has an idea?
(edit: I added a second example)


texts: Average Nearest Neighbor tool extreme restuls for different area?

texts:
I'm using the "Average Nearest Neighbor" tool in ArcGIS 10.3 to calculate whether my point feature class is clustered, random or dispersed.

One of the parameters is area. I tried putting different areas and discovered that if I cut the area lets say by 50% the result changes form clustered (one extreme of the graph created by the software) to dispersed (the other extreme of that very same graph).

Can someone please explain how area is being used by ArcGIS to calculate that? As far as I knoe Nearest Neighbor has nothing to do with area it's all about the distance between the points.


texts: Unable to display geojson data on a leaflet map in a Django app

texts:
I'm developing a Django application using Django-leaflet and Postgresql / PostGIS. I retrieve the data in views.py and I convert them in GeoJSON format  using Django's serializer.
My views.py:

from geonode.geoloc.models import Transport
from django.template import RequestContext
from django.core.serializers import serialize
import json, ast

class LookupView(FormView):
    template_name = 'geoloc/lookupresults.html'
    form_class = LookupForm

    def get(self, request):
        return render_to_response('geoloc/lookup.html',  RequestContext(request))

    def form_valid(self, form):
        # Get data
        latitude = form.cleaned_data['latitude']
        longitude = form.cleaned_data['longitude']

        # Look up roads
        roads =  Transport.objects.all()[0:5]

        roads_json = serialize('geojson', roads,
               fields=('geom',))

        # I remove crs key because it is not recommended by http://geojsonlint.com/
        new_roads_json = json.loads(roads_json)
        new_roads_json.pop('crs', None)
        new_roads_json = json.dumps(new_roads_json)

        # Render the template
        return self.render_to_response({
                          'new_roads_json': new_roads_json
                         })

My template when the form is valid:

{% extends "geoloc/geoloc_base.html" %}
{% block content %}
{% load leaflet_tags %}
{% leaflet_js %}
{% leaflet_css %}

{% endblock %}

And my models.py:

from django.contrib.gis.db import models

class Transport(models.Model):
    gid = models.AutoField(primary_key=True)
    aa = models.SmallIntegerField(blank=True, null=True)
    name = models.CharField(max_length=100, blank=True, null=True)
    geom = models.MultiLineStringField(srid=2100, blank=True, null=True)
    objects = models.GeoManager()

    class Meta:
        managed = False
        db_table = 'transport'

    def __str__(self):
        return self.name

    def __unicode__(self):
        return self.name or u''

When I send the data while geojson content is valid, nothing is displayed on the leaflet map.
Also, no errors in browser's developer tools / firebug or in my django development server.
To test if geojson object that serializer creates is valid, I sended the data as a text with:

document.getElementById("id").innerHTML = geojsonFeature

That geojson text, I posted in http://geojsonlint.com/ and the geometry was displayed on the map with no errors.

Any suggestions ?


texts:
I want to calculate the general orientation of a multilinestring, I used this script to compute the orientation of a linestring :

ST_Azimuth(  ST_Transform(ST_SetSRID(ST_Point(ST_XMax(OSM_geom),ST_YMax(OSM_geom)),2154),432), ST_Transform(ST_SetSRID(ST_Point(ST_XMin(OSM_geom),ST_YMin(OSM_geom)),2154),4326)));

but I don't know how to generate the global orientation of a multilinestring.

Simple Orientation :

Global Orientation :


texts: Automated search of misplaced points with Python or Model Builder

texts:
I am using ArcGIS 10.2.2.

My work is to find points which are misplaced.

I want to know if there's a way to automate this using Python or Model Builder?

I wish to place this in the attribute table which will output a '1' if there's a misplacement of a point in a field.


texts:
Is there a button in QGIS that will allow me to create parallel lines?  I do not mean by using a script, or distance from.

I usually use ArcMap, and in ArcMap you just begin editing, select parallel lines, and click on the line you want to parallel from.  There is also a perpendicular and right angle option. I use these a lot at work to create "perfect" lines around buildings.

Can you do this in QGIS, and if so how do you do it?

The question, "Creating parallel lines at particular distance from lines which are part of a shapefile," does not answer my question. I do not want to do distance from.  In fact I do not want to parallel form line features.  I want to parallel form polygon Features.

Say I am drawing a square with a polygon, I want to be able to make the first line of the square, then right angle from that line, and then right angle until I complete the square.  (Making a perfect square)

Now that I have made a square around a building, I need to draw the boundary of the site.

So I would start by paralleling off of the building lines.

I do this at work on ArcMap all day.  It creates beautiful maps for my clients.  I am trying to find out if I can do this on QGIS.

Again, this is editing polygons, not lines. In ArcMap this editing tool is a toolbar called ****Feature Construction****


texts: How do I show metrics (e.g., population density, counts, literacy rate) based on a custom polygon/shapes?

texts:
I currently have the GeoJSON files containing the census/administrative polygons and metrics associated with them. However, I want to add the ability to view that same data with custom polygons. For example:

I want to drop a marker and view the metrics in a 5km radius
I want to draw a rectangle and view the metrics only in that area

Could you please recommend the best way to approach this problem - right now I am thinking along the lines of finding the % intersection of the custom polygon with the underlying administrative polygon to get the proportionate metrics. Is this even feasible?


texts: Projecting ArcGIS feature dataset containing relationship classes

texts:
I have a geodatabase with more then 50 feature classes and 60-70 tables.
The feature classes are inside a feature dataset which is using a projected cooridnate system.
There are 250+ relationship classes inside the feature dataset linking feature classes with each other as well as with tables.
Similarly there are relationship classes outside feature dataset as well.

My task is to convert the projection of all feature classes from Projected to WGS84.
If I do it one by one and write to a new feature dataset I lose the relationship class.
If I apply project tool on feature data set then it creates a new feature dataset with feature classes but only few relationship classes are there.

How should I project a feature data set (all feature classes)  keeping all the relationship classes?

Even if a break all relationship classes and create a new feature data set is there a way that I can copy the relationship classess without creating each one by one?

I am using ArcGIS Desktop 10.5 and I have ArcGIS Pro 2.3 as well.


texts: Postgis 2 - ST_area : huge difference between 3857 and 4326

texts:
I've been doing tests between projections with Postgis 2.3 and I was surprised by the result of the ST_area function with the Web Mercator projection (3857).

Here is a sample query and result on a PostgreSQL 9.6 database and Postgis 2.3 :

SELECT code,
       name,
       st_area(geom_wgs84::geography))/1000000 AS surface_km2_4326,
       st_area(st_transform(geom_wgs84,2154))/1000000 AS surface_km2_2154,
       st_area(st_transform(geom_wgs84,3857))/1000000 AS surface_km2_3857
FROM communes_www
WHERE numero = '45234'

Here is the result for the city of Orléans, France :

CODE    NAME      surface_km2_4326      surface_km2_2154    surface_km2_3857
------------------------------------------------------------------------------
45234   ORLEANS   27.624689491386977    27.587944000075932  61.377449120370805

As you can see, the area for projections 4326 and 2154 are correct, but the surface for 3857 is completely off.

Is this a bug in Postgis or am I doing something wrong ?


texts: Summing Raster cells to equal circles with specific sum using ArcMap?

texts:
I have application user density data that I've collected. Each gps point has a user density associated with it (for example: 150 users within a 1 km radius, which calculates to users/[pi*R-sq] OR 47.746 users per square km).

I have taken these points and turned them into an interpolated raster surface with a cell size of 50mx50m. Now I can see the whole user density across the map. I accounted for the cell size by dividing the km-sq density by 400 (there are 400 50m by 50m cells in 1 square km). So now each cell reports the number of users (like .25 users per sq 50m).

Now what I want to do is tell ArcMap to create circles that maximize coverage of the map and have within them a total of 150 users. So I want ArcMap to sum up the cell values to equal 150 but within a circular area.
I also want to minimize the number of circles that are needed to cover the whole map.

I'm stuck. I've tried block statistics, but I can't designate that I want users to be added together, only the size of the areas. At then end of this I want a map with circles of 150 users each.

Any ideas?


texts: Qgis freezes when attempting to save geoprocessing user script

texts:
I have recently installed QGIS 2.18.20 32 bits (QGIS-ltr) on Windows 10 64bits machine. I have a small python script in the following path:

C:\PracticaQGisPlugins

The script its visible too in the processing toolbox > scripts > user scripts.

While it executes well, when I first tried to save a modification, an error message showed up:

An error has occurred while executing Python code:

RuntimeError: maximum recursion depth exceeded
Traceback (most recent call last):
    File "C:/OSGeo4W/apps/qgis-ltr/./python/plugins\processing\gui\ProcessingToolbox.py", line 163, in updateProvider
        item.refresh()
   File "C:/OSGeo4W/apps/qgis-ltr/./python/plugins\processing\gui\ProcessingToolbox.py", line 377, in refresh
        Processing.updateAlgsList()
   File "C:/OSGeo4W/apps/qgis-ltr/./python/plugins\processing\core\Processing.py", line 184, in updateAlgsList
        Processing.reloadProvider(p.getName())
   File "C:/OSGeo4W/apps/qgis-ltr/./python/plugins\processing\core\Processing.py", line 189, in reloadProvider
        algList.reloadProvider(providerName)
   File "C:/OSGeo4W/apps/qgis-ltr/./python/plugins\processing\core\alglist.py", line 55, in reloadProvider
        p.loadAlgorithms()
   File "C:/OSGeo4W/apps/qgis-ltr/./python/plugins\processing\core\AlgorithmProvider.py", line 56, in loadAlgorithms
        self._loadAlgorithms()
   File "C:/OSGeo4W/apps/qgis-ltr/./python/plugins\processing\modeler\ModelerOnlyAlgorithmProvider.py", line 53, in _loadAlgorithms
        self.algs = [CalculatorModelerAlgorithm(),
   File "C:/OSGeo4W/apps/qgis-ltr/./python/plugins\processing\core\GeoAlgorithm.py", line 52, in __init__
        self._icon = QIcon(os.path.dirname(__file__) + '/../images/alg.png')

RuntimeError: maximum recursion depth exceeded

Python version: 2.7.4 (default, Apr  6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)]  QGIS version: 2.18.20 Las Palmas, e8e15d5b1e
Python Path: C:/OSGeo4W/apps/qgis-ltr/./python/plugins\processing C:/OSGeo4W/apps/qgis-ltr/./python C:/Users/User/.qgis2/python

Now, the message is not showing up but QGis freezes everytime I try to save the script.

Is this related to the bug mentioned in the following post, which implies I have to reinstall another older version of QGis?:

Configuring QGIS 2.18 processing script locations

What can be done?


texts: From CSV file (latitude and longitude fields) to point .shp

texts:
How i create a new point shapefile from CSV table?

If i have a CSV file with a latitude and longitude fields, is there a method for to have in output a puntal shapefile?


texts: Unexpected behavior of collect() in geometry generator

texts:
I'm having unexpected results using geometry generator with the following expression:

 buffer(collect($geometry),0)

…which was an answer given in this topic. I have used this previously with the expected result, but now something is strange. It appears that the geometry is being collected (and buffered) not once, but over and over—see the attached screencast.
It's not infinite; it will do this about 20-25 times and then stop.

I can't tell which part of the function is the problem. If I just call collect($geometry), I'll get the same result (but rendered faster as it's not buffering). If I just call buffer($geometry, 0), then it renders some features only once, and others a varying multiple of times (see screenshot).

I'm sure it's data related because it doesn't happen on all layers. What can I look for that would cause this multiple rendering?


texts: Change separator of attributions in Leaflet?

texts:
Is there a way to change the comma shown in the image below to |

I used the QGIS plugin QGIS2WEb to create the Leaflet map and wanted to link to the source of the information.

Currently, in the HTML code I have used the attribution format like so:

var overlay_GoogleHybrid_0 = L.tileLayer('https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}', {
        opacity: 1.0,
        attribution:'Base Layer Map info &copy; Google',
    });

I am also unsure if this is the correct method of referencing the Google Hybrid layer I am using (Obtained using the Quick Map Services plugin).


texts:
I'm trying to read multiple .asc files in R so I can perform some analyses later on them. I'm using this code:
library(raster)
list_file <- list.files ("C:\\Users\\Areej\\Documents\\Gibbon_bangladesh\\Bangladesh_Env_Vs_30s_ascii_R\\",
 pattern ".asc", full.names = T)
gibbon <- stack(file_list)
gibbon

but it shows me unexpected token for the ".asc" part in the code


texts:
I'm trying to georeference this part of an old map of Indochina.

.

The full extent gives more information:

.

To do so, I used QGIS and tried to infere the datum, I tried "Unknown datum based on Everest", but it gave me a big difference that you can see in this screenshot:

.

I did it with the three geodesic intersections available in the map.
To understand the problem, I looked at Luang Prabang (the city close to the intersection of geodesics on the top of the map) coordinates and found 19.53 N and 102.80 W on geonames.org. On my map, it is more about 100 W, that is a big difference.

Thus my questions:

given that it is an old map, does the inaccuracy result from the mappers who did the map (hats off for there impressive work btw !) or is it that I do not have the right method (datum, choice of points, ...) ?
how can I georeference this map ? (My goal is to digitize the borders).


texts:
I hope I can get same help from here.
I started to put together a leaflet map, please see here: http://gis.xyz/leaflet2.html
I want to show my postgis point (uk based) features on this map with cluster method. The points are already on the map (I used geoserver- jsonp), but I cannot figure out how can I add to the layer to the cluster layer.
I also added a markercluster example to the script, this is how I would like to see my layer.
It would be great if someone would help me to answer this question.

     var owsrootUrl = 'http://217.8.255.188:8080/geoserver/opengeo/ows';

 var defaultParameters = {
     service : 'WFS',
     version : '2.0',
     request : 'GetFeature',
     typeName : 'opengeo:evernote_geom',
     outputFormat : 'text/javascript',
     format_options : 'callback:getJson',
     SrsName : 'EPSG:4326'
};

var parameters = L.Util.extend(defaultParameters);
var URL = owsrootUrl + L.Util.getParamString(parameters);

var WFSLayer = null;
var ajax = $.ajax({
    url : URL,
    dataType : 'jsonp',
    jsonpCallback : 'getJson',
    success : function (response) {
       WFSLayer = L.geoJson(response, {
            style: function (feature) {
                return {
                    stroke: false,
                    fillColor: 'FFFFFF',
                    fillOpacity: 0
                };
            },
            onEachFeature: function (feature, layer) {
                popupOptions = {maxWidth: 600};
                layer.bindPopup('
'+feature.properties.url+'

'+feature.properties.title
                    ,popupOptions);
            }
        }).addTo(map);
    }
});


texts: Setting Geowebcache to a long time expiration

texts:
I have a Disk blobstore to cache my WMS layers. Fine.
Now I need to set it to a very long ( a month? ) time to expire the disk cache because I have just a few users using and is not enough to keep my tiles cached.

How can I do this?


texts: Cost distance of one raster using routes of separate cost distance in QGIS

texts:
Is there a way, using QGIS, I can calculate the accumulated for one raster of cost values but instead of using the least cost distances from that raster, using the routes/paths determined from a separate cost raster?
For instance, say I have a cost raster for energy expenditure based on slope but I want energy costs taking the least cost distances for another raster such as the vegetation density. The product would essentially be the sum of said energy values but while walking through the thinner vegetation if possible.
I need this for an entire surface rather than a path.


texts: I have black pixels inside water and land. I have covered the land using a shapefile. How do I mask out the entire land so that no pixel is counted on the land? I only want to count pixels inside the water.


texts: Combining 2 rasters using Con and IsNull in Python

texts:
I am trying to write a Python script that converts most NoData pixels to 0 in a raster (Raster1 below), except for pixels that overlap with a pixel value of 1 within another raster (Raster2, which has pixels that are either 0 or 1). I am therefore trying to combine Con and IsNull in the following statement, where Raster1 is the raster whose NoData pixels should be converted to 0, and Raster2 is raster containing the values 0 and 1:

Con(((IsNull(Raster1)) & (Raster2==0)), 0, Raster1)

The output raster is correctly converting NoData pixels to 0, but it’s not completing the 2nd part of the condition, which specifies that a pixel in raster1 should only be converted to 0 if both the raster1 pixel is NoData and Raster2 pixel=1.

I have a workaround where I develop the Con statement without the Raster2 portion of the argument, and then multiply this output by a mask. However, I worry that this would take a lot longer to process, as the rasters I’m working with are huge. I’m also trying to understand how to use the Con statement more, and I’m confused why this is not working.

My rasters are TIFFs, although I do perform a step earlier in the script that converts a mosaic to a TIFF. I don’t see how this would be causing the problem, but I’ve included the full script below just in case (the Con statement is in the last section of the script):

import arcpy, os
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")

# Variables
inWS = r"G:\ellie\_testing_script\mosaics.gdb" # UPDATE
outWS = r"G:\ellie\_testing_script\out_TIFs" # UPDATE
Raster2 = r"G:\ellie\MosaicExports\GAP\GAPStatus12_All" # UPDATE
scratch = r"G:\ellie\_testing_script\scratch" # UPDATE
print ("Variables Set")

# Environments
env.workspace = inWS
arcpy.env.snapRaster = "G:\\ellie\\SnapRaster\\NLCD_2011_USGS_LandCover_AlbConEA_16bitU.tif"
env.overwriteOutput = True print ("Environments Set")

# Create list of mosaics in Geodatabase
mosaics = `arcpy.ListDatasets("*","Mosaic")`

# Loop through mosaics in geodatabase: for mosaic in mosaics:
    print("Processing " + mosaic)
    # We aren't using the GAP123 inputs, so don't process those:
    if mosaic.startswith("RSR_GAP123"):
        pass
    else:
        # Loop through mosaics (not GAP123) to export to new TIF raster
        inMosaic = inWS + "\\"  + mosaic
        # outName needs to be determined by the 2nd "_" argument in order to change some output filenames
        second_argument = mosaic.split('_')[1]
        last_arguments = mosaic[10:]
        if second_argument == "GAP12":
            outName = "PWRSR_GAP12_" + last_arguments + ".tif"
        else:
            outName = mosaic + ".tif"
        arcpy.MosaicToNewRaster_management(inMosaic, scratch, outName,"", "32_BIT_FLOAT", "", "1", "LAST", "FIRST")
        print(outName + " exported as TIF. Converting NAs to 0s...")

        # Convert NA to 0s, except for PAs
        Raster1 = scratch + "\\" + outName
        outCon = Con(((IsNull(Raster1)) & (Raster2==0)), 0, Raster1)
        print(outName + " NAs converted to 0s. Saving raster...")
        outCon.save(outWS + "\\" + outName)

print ("Finished")


texts: How to clip a RGBA raster with another single band raster using QGIS or GDAL

texts: I have a RGBA raster and I would like to set all the data over the continents to transparent. I have a single band raster where land is equal to 0 and ocean is equal to 1. Normally if I had a single band raster, I would multiply the single band raster with the above mask raster, in raster calculator, and this would remove the land data. How do I do the same thing with a RGBA raster in QGIS or GDAL, where the result would be 0, 0, 0, 0 for the four bands?


texts: QGIS Opacity issue with polygons and labels

texts:
I have a layer with a number of polygons. I can happily adjust the opacity of the polygons using the normal slider, but when I add labels to the polygons from the Attributes table the labels do not change their opacity in sync with the polygons.

Have I missed a setting or does this need an expression/code somewhere to create the linkage?


texts: Converting .kmz/.kml to .shp using QGIS?

texts:
I have a file of 3d buildings created in SketchUp that I need to convert to shape format. Tried to export the file as .kmz and importing it to QGIS for the conversion, but when the .kmz file opens in QGIS, it contains no features.

The same happens if I unzip the kmz to kml before opening it in QGIS. I have tried to figure out if there is a plugin that will help me but without success. Im working on QGIS 2.18.11.

My original file is too large (<70 mb) to use any of the online conversion sites.

I have previously imported .kml files in QGIS without any problems, so I'm suspecting the error is in the files that I export from SketchUp.

Does anyone know what could be causing this or have any suggestions of how to proceed?


texts:
Is there a way to Order the way Layers are displayed when adding them to a process in the modeler?

I am trying to add 4 layers to a VRT process before using a Clip by polygon function and the problem I am having is this-

When adding the layers to the VRT process they are coming in as

Raster 1,
Raster 3 (optional),
Raster 4 (optional),
Raster2.

This is annoying as I would like them ordered 1,2,3,4 as layers 3  and 4 are optional and wont be used each time when the model is being used. I just can't find how to set which order they are asked for in the execution of the model.


texts: What are the pros/cons of using an MSSQL SDE Database rather than a File-gdb in ArcGIS Server

texts:
I have taken over an ArcSDE 9.3 Enterprise database setup which has around 80% of it´s vector data in an MSSQL SDE Database system and 20% in FileGDB on the application server. I am trying to understand if there is a logical reason as to why someone would set this up. Can anyone offer reasons why a MSSQL SDE Database is prefered to FileGeodatabases?  Versioning is not being used in either case. The web clients which have read/write access seem to all be in the MSSQL SDE DB. There are also 6 SDE Services running although most of the connects are direct.

If anyone could explaing the logics behind this (if there are any) I would be grateful.


texts: How to calculate Median of a table field in the Dissolve Tool of ArcGIS?

texts:
The Dissolve Tool of ArcGIS 10.0 has some statistical parameters as a type of calculation  new values of each field in a process of dissolving of identical parameters . It has SUM, MAX, MIN, RANGE, STD and a bit more Statistical types , but it doesn’t have MEDIAN.

It’s very  strange, why it doesn’t have  Median (and some others important statistical parameters). But maybe it’s possible to “turn on” calculation  “Median” in the Dissolve Tool, or solve this task in some other way?