texts:
This resource shows how to use NULL values in pyqgis. Within the python colsole the comparison value == NULL passes without exception. However, when run from my plugin NULL is neither in the global nor the local namespace

NameError: global name 'NULL' is not defined

The only way to check for NULL values seems to be to check for an instance of QPyNullVariant

from PyQt4.QtCore import QPyNullVariant
if isinstance(value, QPyNullVariant):
    pass

How do I import NULL so I can compare value == NULL?

EDIT

For those wondering why and how NULL is different from None have a look at this blog entry. Essentially historically is wasn't possible to return None:

  Turns out by removing QVariant from PyQt it had some impact on methods that expected a NULL QVariant - A QVariant with no value. Passing None didn't work because those methods needed the type information that QVariant holds, even when NULL.

On None comparisons it says:

  One way to check if something is Null in Python is to use value is None however this will not work with our NULL type. Overloading the is operator in Python is not supported and there is no way we can support this - trust me I have tried. is is really doing id(object) == id(object) under the hood:


texts: When will QGIS use the WMS supplied legend images ? is there a workaround?

texts:
I use a lot of wms services with complicated maps : a legend is essential.
Qgis displays none of those legends. The URL is quoted in the metadata.

example
WMS

http://geo.agiv.be/ogc/wms/product/ANB?/Grnkrt13

layer "Groenkaart Vlaanderen 2013"

legend url :
http://geo.agiv.be/ogc/wms/product/ANB?request=GetLegendGraphic%26version=1.3.0%26format=image/png%26layer=Grnkrt13


texts: How to set margin on PageLayout in ArcGIS Engine?

texts:
When we create a PageLayout from the A4 Landscape template and print the layout on a pinter with A4 Landscape paper we get four pages instead of one. We are using the IPage.PageToPrinterMapping = esriPageMappingTile and the reason that we get four instead of one page is that the printer by default has a margin for the A4 format. If I try to do the same procedure for a format that the printer doesn't have a margin for it prints only one page. Is there a way to create the same margin on the PageLayout?

Steps to reproduce in ArcMap:

File | New...
Choose Templates | ISO A4 Landscape
File | Page and Print setup...
Set Printer and Map to the same
format A4 and orientation landscape
File | Print...
The Print dialog is showing four pages

Update: Tried different approaches using the IPage interface but can't find any API to set the margin.

Tried setting the PrintableBounds on the page, but it seems to be readonly:

m_hookHelper.PageLayout.Page.PrintableBounds.Envelope.PutCoords(docPaper.PrintableBounds.XMin, docPaper.PrintableBounds.YMin, docPaper.PrintableBounds.XMax, docPaper.PrintableBounds.YMax)

Tried setting the Custom size on the Page and it sets the size correctly, but it corp  the layout, because it doesen't set the offset:

m_hookHelper.PageLayout.Page.PutCustomSize(docPaper.PrintableBounds.Width, docPaper.PrintableBounds.Height)

Tried to set the FormID to be esriPageFormSameAsPrinter and it workes as long as the user chooses the layouts' orientation and format on the printer. If not, it crops the layout when it is printed.

m_hookHelper.PageLayout.Page.FormID = esriPageFormID.esriPageFormSameAsPrinter


texts: How to check if a point is within bounded rectangle from Java

texts:
I make an AJAX request to retrieve data from a mysql database. I add the data to google.maps.Data layer as geojson. I would only like to load the data based on the viewport. To obtain the viewport bounds I use:

var bounds = map.getBounds();

With the bounds information, I have sent an AJAX request to my server to get the markers that are within those bounds.

I query the database table for all the data it has, can I instead query mysqlonly for the data contained within the bounds. If that is not an option, how do I check if a point is within the rectangle (on the server)?


texts: Google Earth Engine Failed to decode JSON

texts:
var S2Collection = ee.ImageCollection('COPERNICUS/S2')
    .filterBounds(geometry)
    .filter(ee.Filter.calendarRange(10,02,'month'));

var cloudBitMask = ee.Number(2).pow(10).int();
var cirrusBitMask = ee.Number(2).pow(11).int();

function MaskS2Clouds(image) {
  var qa = image.select('QA60');
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
             qa.bitwiseAnd(cirrusBitMask).eq(0));
  return image.updateMask(mask);
}

var S2CloudMasked = S2Collection.map(MaskS2Clouds);

function addNDVI(image) {
  return image.addBands(image.normalizedDifference(['B8', 'B4']).rename('NDVI'));
}

var unfNDVI = S2CloudMasked.map(addNDVI).select(['NDVI']);

function filterNDVI(image){
  var masked = image.updateMask(image.gt(0.2));
  var mean = ee.Number(masked.reduceRegion({
      geometry: geometry,
      reducer: ee.Reducer.mean(),
  }).get('NDVI'));
  return image.addBands(masked.divide(mean.getInfo()).rename('stdNDVI'));
}

var fNDVI = unfNDVI.map(filterNDVI).select(['stdNDVI']);

My code fails on the last function. I want to divide each masked image by its mean (pixel by pixel) but I get:

  Failed to decode JSON.
  Error: Field 'value' of object '{"type":"ArgumentRef","value":null}' is missing or null.
  Object: {"type":"ArgumentRef","value":null}.

What is wrong?


texts:
Is it possible to import styles in ArcMap?

Background:

We made a workflow in FME where CAD files are written to shapes and sld's which can be loaded in QGIS. This way we can view CAD data in QGIS as it was drawn in CAD. This really bridges the CAD - GIS gap.

Now we want to deploy this data AND styling to ArcGIS Online. After some digging around I think I understand how styling in ArcMap works:

Define style in style library.
Link style to featureclass.

Because it is designed this way it might be hard to import styles like I do in QGIS. I wonder if someone has worked around this challenge? We're talking about a serious bunch of cad files from all types of file owners. So no default styling / layers / etc.

The best alternatives I can think of:

WMS QGIS Server for style combined with featureservice for vector.
Generating tiles with QGIS QTile and upload to AGOL.
Generating raster with FME and tile with ArcMap.

But I would like to only use a featureservice.

Update 201610006:
Tiles are out of the question because of processing time and storage.

The logic step is to map the featureclasses to company standards with standard styles but that is not bridging the CAD-GIS gap.

Brainstorming about this we see another possibility: Writing to service definition file (zip with json) and edit this one before uploading. SDF contains data AND styling. I think it might be a long shot and deep dive in json but with FME everything is possible. The downside is the possibility ESRI closing this door after reading my posts about it.


texts:
I need to colorized geometries from shapefiles and convert them to raster image (geotiff).

I tried to use geotools's VectorToRasterProcess and ImageWorker
("d:\gis\tmp\32 27 25 w4.shp" is path to my test shapefile):

URL url =  new File("d:\\gis\\tmp\\32 27 25 w4.shp").toURI().toURL();
ShapefileDataStore ds = new ShapefileDataStore(url);
FeatureSource fs = ds.getFeatureSource("32%2027%2025%20w4");

GridCoverage2D raster = VectorToRasterProcess.process(fs.getFeatures(),
     "ID", new Dimension(400, 400), fs.getBounds(),"test", null);

ImageWorker imageWorker = new ImageWorker(raster.getRenderedImage());
imageWorker = imageWorker.rescaleToBytes();
imageWorker.makeColorTransparent(Color.BLACK);
imageWorker.forceComponentColorModel();

raster = new GridCoverageFactory().create("one", imageWorker.getRenderedImage(),  fs.getBounds());

GeoTiffWriter writer = new GeoTiffWriter(output);
GeoTiffWriteParams params = new GeoTiffWriteParams();
ParameterValue<GeoToolsWriteParams> value = GeoTiffFormat.GEOTOOLS_WRITE_PARAMS.createValue();
value.setValue(params);
writer.write(raster, new GeneralParameterValue[]{value});
writer.dispose();

The result geotiff is here

What is not right in my configuration of ImageWorker?
Give, please, some advices.

Regards,
Dmitry


texts: Digital Globe WMTS won't work in MapInfo 12.0

texts:
Digital Globe recently provide API and WMTS service to the public.

I follow the instruction in the https://mapsapi.digitalglobe.com/ and get the token here https://mapsapi.digitalglobe.com/account

In QGIS adding WMTS work like this: open the WMS dialogue and put the API link together with token ID provided by Digital Globe.

  http://api.mapbox.com/styles/v1/digitalglobe/ciode6t5k0081aqm7k06dod4v/wmts?access_token="token
  here"

In QGIS it work fine.  But in MapInfo when I put the same URL above to MapInfo WMS dialogue I get the error "Http status 422. Enable to retrieve capability form server".

Any one has successfully adding Digital Globe WMTS to MapInfo?


texts: Where I can find ArcGIS online symbols to use in Desktop

texts:
I'd like to use this symbols:

In desktop ArcGIS.

Is this possible?


texts: Using BOTH size and colour to represent 2 pieces of information on single point in QGIS

texts:
I want the colour of the point to represent value from data column x (eg. rating from 1-5) and SIMULTANEOUSLY the size of the point to represent value from data column y (eg. Number of ratings 0-600).

I can swear that I have done this before but now I can't find how to do this in QGIS 2.10. It only gives the options to display using colour OR size.


texts: How to remove ghost lines around polygons in QGIS 2.12?

texts:
I'm trying to remove the very thin, yet visible white "ghost" lines in between polygons. You can see This in the image below.

Hopefully you don't have to strain your eyes to see that there are very thin white lines in between each of the polygons. I've made the lines transparent and given them a thickness of 0, but they're still there. Here are my settings:

Any thoughts on how to solve the problem or workaround solutions?


texts: JPEG file wont open in QGIS 3.10 but does in QGIS 3.4

texts:
I recently made the upgrade to QIGS 3.10 from 3.4 and I am having trouble loading a JPEG image that previously worked fine in 3.4. In the georeference window when I try to open the raster I am given the following error:
D:/DJI_0151.JPG is not a supported raster data source.
libjpeg: Bogus marker length
If I try to open it as a raster in the map I get the error 'unable to open GDAL dataset'
It is funny because some JPEGS seem to work fine and some don't now when previously they all worked fine.
Any idea what is causing this?


texts: Mapbox tiles not showing on Leaflet map using Google Chrome

texts:
https://underjollyroger.github.io/EarlyIllinoisLandcover
There is a 403 (Forbidden) error and the Mapbox tiles will not display when viewing from Google Chrome (zoom in six times with the mouse wheel to view tiles). I'm using a Mapbox API access token and the app is hosted from Github. The app works fine from Firefox.


texts: Mosaicking layerstacked images in Google Earth Engine

texts:
I am trying to layerstack two set of mosaicked images, or I have two sets of individual layerstacks of seperate swaths which i want to mosaic.

I am attaching the code here:

var roi = ee.Geometry.Rectangle(68.8602, 22.6003,70.0577, 21.6079);

var collection1 = ee.ImageCollection(['COPERNICUS/S2_SR/20200205T055029_20200205T055653_T42QWK',
                                      'COPERNICUS/S2_SR/20200116T055159_20200116T060121_T42QWK',
                                      ]).select(['B8','B4','B3','B2']).filterBounds(roi);
print('Number of images in collection1:', collection1);

var stackCollection1 = function(collection1) {

  var first1 = ee.Image(collection1.first()).select([]);

  var appendBands1 = function(image, previous) {
    return ee.Image(previous).addBands(image);
  };
  return ee.Image(collection1.iterate(appendBands1, first1));
};
var stacked1 = stackCollection1(collection1).clip(roi);
print('stacked image1', stacked1);
Map.addLayer(stacked1, {min:0, max:4000}, 'stacked1');

var collection2 = ee.ImageCollection(['COPERNICUS/S2_SR/20200205T055029_20200205T055653_T42QVK',
                                      'COPERNICUS/S2_SR/20200116T055159_20200116T060121_T42QVK'
                                      ]).select(['B8','B4','B3','B2']).filterBounds(roi);
print('Number of images in collection2:', collection2);

var stackCollection2 = function(collection2) {

  var first2 = ee.Image(collection2.first()).select([]);

  var appendBands2 = function(image, previous) {
    return ee.Image(previous).addBands(image);
  };
  return ee.Image(collection2.iterate(appendBands2, first2));
};
var stacked2 = stackCollection2(collection2).clip(roi);
print('stacked image2', stacked2);
Map.addLayer(stacked2, {min:0, max:4000}, 'stacked2');


texts: pgr_astar function does not accept negative values?

texts:
I'm using postgreSQL and I'm running the pgr_astar function to get the optimal path. I calculated my path using a formula that uses angle of elevation or the angle of depression as input. So I'm getting negative values on the cost column, I pgr_astr function does not return any result when I'm using the cost values.

Any suggestion please? How can I manipulate this function to get the optimal path using the cost's negative values?


texts:
I have three tables containing increasingly-complex geometry. I'm clipping Tables 2 and 3 to the first, least-complex table (with 700 multipolygons). The other two tables contain spatial attributes I'm after (t2: 300 000 records, t3: millions).

Table 3 - most complex - has less than 100% areal coverage, so I need to "see through" the holes to the geometries of Table 2.

With SO help sorted a couple intersect queries, read docs explaining indexing, vacuum / explain, and came back to speed these up:

SELECT ST_Intersection(t1.geom, t2.geom)
FROM t1
JOIN t2 ON ST_Intersects(t1.geom, t2.geom)

... when I found this thread from 6 years back, a massive help, taking initial queries from 30 minutes to under 4:
Acquiring ArcGIS-like speed in Postgis

Using the ST_CoveredBy function and CASE flow we have two nearly-identical queries relating t3 with t1, then t2 with t1:

SELECT t1.location, t3.attributecode,
  CASE
  WHEN ST_CoveredBy(t3.geom, t1.geom)
  THEN t3.geom
ELSE
  ST_Multi(
      ST_Intersection(t3.geom, t1.geom)
      ) END AS geom
FROM table3 AS t3
    INNER JOIN table1 AS t1
    ON ST_Intersects(t3.geom, t1.geom)
    );

I've made tables from these, here's t2 cut to t1 boundaries:

With t3 overlying this (transparency added to make it clear they're two geom tables):

Next is aggregating t2 and t3, using ST_Union perhaps, which is slow even on smaller test datasets.

What are some ways of speeding this up?


texts: How to add FP_WGS_84_IBCAO_Polar_Stereographic projection to QGIS?

texts:
I have a GeoTIFF image that is projected (see below), I would like to add this projection to the list of CRS in QGIS.

Can someone help me do this?

PROJCS["FP_WGS_84_IBCAO_Polar_Stereographic",
   GEOGCS["FG_WGS_84",
      DATUM["FD_WGS_1984",
         SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],
      PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
      UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],
   UNIT["meter",1,AUTHORITY["EPSG","9001"]],
   PROJECTION["Polar_Stereographic"],
   PARAMETER["latitude_of_origin",75],
   PARAMETER["central_meridian",0],
   PARAMETER["scale_factor",1],
   PARAMETER["false_easting",0],
   PARAMETER["false_northing",0],
   AUTHORITY["EPSG","3996"],
   AXIS["X",
      UNKNOWN],
   AXIS["Y",
      UNKNOWN]]


texts:
I've come across the following coordinates:

x = 600907.76
y = 2437133

It's a coordinate located in Hawaii. I thought it was UTM, but the negative x coordinate makes me think otherwise. Anyone know what it should be?


texts: Creating geodata service from version other than DEFAULT in ArcGIS for Server?

texts:
This is with ArcGIS for Server and Desktop 10.1

I want to host a geodata service, but I don't want it to be from the DEFAULT version. I want it to be from an alternate version, that I can check the changes on before I post it to the DEFAULT.

I don't see any options to do this.

Is it possible?


texts: Using QCompleter with strings composed by different words in QGIS Plugin - PyQT5

texts:
I'm using QCompleter in my QGIS plugin. The code below works well, but I'd love to make it more sophisticated. In this example I have some strings that are composed of some words, for example "How to use QCompleter". In this moment, to see the completer I need to start writing by the first word ("How"). I'm wondering if there's a way to see the completer also writing in the LineEdit any word in the string, for example "to" or "use".

strList = ["Prove PyQt5", "How to use Qcompleter", "I love StackExchange"]
completer = QCompleter()
completer.setCaseSensitivity(0)
dlg.lineEdit.setCompleter(completer)
model = QStringListModel()
model.setStringList(strList)
completer.setModel(model)


texts: Loading shapefile in GeoDjango and displaying it in admin section?

texts:
I am trying to load a .shp file in GeoDjango Admin. As a newbie I followed their official documentation and based on that I am trying to load my own shapefile. No errors are being shown but I can't see my shapefile being loaded. All I can see is a world map but I can see attributes related to my shapefile.
Here is the code:
load.py

from django.contrib.gis.utils import LayerMapping
from .models import TestGeo

testgeo_mapping = {
    'state_n': 'STATE_N',
    'admin_stat': 'ADMIN_STAT',
    'geom': 'MULTIPOLYGON',

}

state_shp = os.path.abspath(
    os.path.join(os.path.dirname(__file__), 'data', 'STATE.shp'),
)

def run(verbose=True):
    lm = LayerMapping(TestGeo, state_shp, testgeo_mapping, transform=False)
    lm.save(strict=True, verbose=verbose)

models.py

class TestGeo(models.Model):
    state_n = models.CharField(max_length=50)
    admin_stat = models.IntegerField()
    geom = models.MultiPolygonField(srid=4326)

    def __str__(self):
        return self.state_n

Output:

Output of python manage.py ogrinspect forest/data/STATE.shp TestGeo --srid=4326 --mapping --multi:

class TestGeo(models.Model):
    state_n = models.CharField(max_length=50)
    admin_stat = models.IntegerField()
    geom = models.MultiPolygonField(srid=4326)

# Auto-generated `LayerMapping` dictionary for TestGeo model
testgeo_mapping = {
    'state_n': 'STATE_N',
    'admin_stat': 'ADMIN_STAT',
    'geom': 'MULTIPOLYGON',
}

What am I doing wrong?

shapefile is completely fine I have plotted it using matplotlib and also inspected it in QGIS.


texts: Japanese Coordinates Not Aligned Correctly

texts:
I am attempting to plot coordinates I obtained on the Nara National Research Institute for Cultural Properties, or Nabunken, website in QGIS. They have a archaeological site search page (http://mokuren.nabunken.go.jp/Iseki/index.html), and they list coordinates for the sites. After putting together a group of sites, I noticed that their positions were all off. I tried inputting the coordinates in Google Earth as well and had the same issue. They all seem to be off by a similar amount, about 30-35km south and 20-25km east of where they should be.

For an easy example the coordinates listed for the Sannai Maruyama site are 40.4845, 140.4158, but bringing up the site on Google Earth gives me coordinates of 40.811062, 140.695648. The website only lists that the coordinates are WGS.

Any ideas of what might be going on?

Here is a copy of the data sheet I am getting the point coordinates from.

The coordinates are circled and the WGS designation title is underlined. One additional bit of possibly useful information is that there are no decimals included for the Sannai Maruyama site coordinates, so the original data shows 404845, 1404158 and on other sites I have only seen decimals before the last 2 numbers, so for another specific section of the site they have the location listed as 404849.6, 1404207.3 .


texts: Looking for ready tools to display real-time positions online

texts:
we want to show that in a big city like Palermo, cycling can be the fastest tool (and at the same time the most ecological). We will give 4 smartphones to 4 users, and there will  be 4 ways to move: by bicycle, by public transport, by foot and by car. We would like to use some smartphone apps to show in real-time in a webpage, in a single map, the route of each user.
I know I can use something like Google Latitude or Glympse to track in real-time user position, but I do not know how to put in a single map 4 real-time different routes.

Our goal is to promote the use of bicycle, because we have a high air and noise pollution and the motor vehicle traffic makes mobility very problematic. We will do it voluntarily without any economic budget, then we are looking for some simple mush-up solution that does not force us to produce several hours of development. Something almost ready for "production".

Do you know some solution?


texts: How would I determine my search radius when conducting a Topographic Position Index in SAGA GIS?

texts:
When running the TPI tool in SAGA, I get this parameter-entry window:

Under the section highlighted in the red box, what do these parameters mean? Do the search/neighbourhood grouping parameters mean cells? What would the "standardize" section indicate? Would distance weighting also be recommended?


texts: Query building - combining multiple fields and values

texts:
I'm using QGIS 2.18.16 and i'm learning to use the query builder to create my own layers. I have noticed two separate fields (name, building) which have similar category of values such as "schools, train station, shops, libraries". I'm seeking to have one query that combines both fields, allowing the creation of one layer which I can then categorise.

How do I create this query?


texts: Automatic rotation of maps in composer to maximize scale of object of interest

texts:
We are working on a composer atlas of hundreds of maps showing highway routes, all of which are north-oriented. We would like to be able to rotate these maps to the ideal rotation that maximizes the scale that the route can be displayed. I believe the best way to calculate this would be by looking at the dimensions of the minimum bounding box of the route, compared to the dimensions of the map frame.

For example, suppose that I have a route whose MBB is 300m in the X dimension and 100m in the Y dimension (3:1 ratio). In comparison, the ratio of the map frame is 2:3 (taller than wider). Is there a way to calculate the angle that the map needs to be rotated away from north so that the MBB of the route at this rotation becomes 2:3, which should maximize the scale that it can be displayed?


texts: Set default Map with saturation in Earth Engine

texts:
When setting the map options to:

Map.setOptions( 'HYBRID', {'HYBRID':[ { featureType:'all', stylers:[{saturation:-100}] },]  }  );

It creates the map as an option to be clicked in the user interface but does not set it as the default map.

How do I set it to be the default display map style?


texts: Calculate area of a raster for an irregular polygon in PostGIS

texts:
I have a number of polygons with irregular shapes (i.e. non-rectangular) and want to calculate their areas after I have converted them to rasters.

Here are my steps:

Convert my geometry to a raster, using a reference raster with fixed pixel sizes (50m x 50m).
Clip the new raster to the old geometry to get rid of the extra pixels outside of the area of interest.
Calculate the area by counting the pixels and multiplying by the known fixed sizes. (This is the part I am having trouble with)

My PostGIS query:

SELECT
    ST_ValueCount(ST_CLIP(ST_AsRaster(my_geometry, reference_raster)), 1) as pixel_count,
    ST_Area(my_geometry) as geom_area
FROM ra

This gives me

 | pixel_count | geom_area
-+-------------------------
0|     (1, 123)|      4567

which is not an array so I cannot index to get the count and directly compare to the geometry area.

I could extract it in Python and but it seems to me there should be a clean way to do this in PostGIS. Does anyone know how to do this?


texts:
What am I doing wrong? I am not able to receive the raster layer. However, when I am opening the url from the browser I am able to receive the image.

from qgis.core import *

url='url=http://geodata.gov.gr/geoserver/ows?service=WMS&version=1.3.0&request=GetMap&bbox=448789,4177572,501682,4222482&crs=EPSG:2100&width=431&height=365&layers=geodata.gov.gr:18fe8202-e9de-4b83-a6ac-d4978b690b46&styles&format=image/png&dpi=96&map_resolution=96&format_options=dpi:96'
rlayer = QgsRasterLayer(url, 'test', 'wms')

print(rlayer.isValid())
# Gives True
QgsProject.instance().addMapLayer(rlayer)

But the layer doesn't load

The warning in the logs is:

  2020-03-11T16:15:41     WARNING    Map request failed [error: Error transferring http://geodata.gov.gr:80/geoserver/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=448789,4177572,501682,4222482&CRS=EPSG:2100&WIDTH=357&HEIGHT=303&LAYERS=geodata.gov.gr:18fe8202-e9de-4b83-a6ac-d4978b690b46&STYLES=&FORMAT=image/png&DPI=96&MAP_RESOLUTION=96&FORMAT_OPTIONS=dpi:96&TRANSPARENT=TRUE - server replied: Service Unavailable url: http://geodata.gov.gr:80/geoserver/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=448789,4177572,501682,4222482&CRS=EPSG:2100&WIDTH=357&HEIGHT=303&LAYERS=geodata.gov.gr:18fe8202-e9de-4b83-a6ac-d4978b690b46&STYLES=&FORMAT=image/png&DPI=96&MAP_RESOLUTION=96&FORMAT_OPTIONS=dpi:96&TRANSPARENT=TRUE]

But the service request in a browser works.

Using the QGIS Network Logger plugin, I can see that after rlayer = QgsRasterLayer(url, 'test', 'wms') A GetCapabilities request is sent:

http://geodata.gov.gr/geoserver/ows?service=WMS&SERVICE=WMS&REQUEST=GetCapabilities' -H 'User-Agent: Mozilla/5.0 QGIS/31002'  --compressed

and a successful response is received.

I guess that the True response from print(rlayer.isValid())
is a reflection of the successful GetCapabilities request.


texts:
I have been asked to export a DEM into DWG format, is this possible? The only thing I can think of is to create a regular grid of points and to export this.

Also I have been asked to save geotiff raster images to DWG is this possible?


texts: Unable to add SHP from GEE Assets into QGIS via QGIS-GEE plugin

texts:
This might come across as a duplicate question here but I cannot find a way out.
I've got two shapefiles uploaded in Assets in GEE against my username.
I am trying to add them into QGIS and it fails all the time.
I can confirm that the GEE plugin has been installed correctly because I can load imagery without any issues. Also, I am able to perform analysis in GEE but want to do so in QGIS because there are certain things that are a bit easier there.
The code I'm trying is as follows:
poly_layer = ee.FeatureCollection('users/myusername/folder/polygon');
Map.centerObject(poly_layer);

The error I'm getting is as follows:
Traceback (most recent call last):   File "C:\Users\myname\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\ee_plugin\extlibs_windows\ee\data.py", line 338, in _execute_cloud_call
    return call.execute(num_retries=num_retries)   File "C:\Users\myname\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\ee_plugin\extlibs_windows\googleapiclient\_helpers.py", line 134, in positional_wrapper
    return wrapped(*args, **kwargs)   File "C:\Users\myname\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\ee_plugin\extlibs_windows\googleapiclient\http.py", line 907, in execute
    raise HttpError(resp, content, uri=self.uri) googleapiclient.errors.HttpError: <HttpError 400 when requesting https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/value:compute?prettyPrint=false&alt=json returned "Collection.loadTable: Collection asset 'users/myname/master_shp/polygon' not found.">

During handling of the above exception, another exception occurred:

Traceback (most recent call last):   File "C:\PROGRA~1\QGIS3~1.14\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)   File "<input>", line 1, in <module>   File "<string>", line 9, in <module>   File "C:/Users/myname/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\ee_plugin\Map.py", line 70, in centerObject
    coords = rect.bounds().getInfo()['coordinates'][0]   File "C:\Users\myname\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\ee_plugin\extlibs_windows\ee\computedobject.py", line 95, in getInfo
    return data.computeValue(self)   File "C:\Users\myname\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\ee_plugin\extlibs_windows\ee\data.py", line 703, in computeValue
    prettyPrint=False))['result']   File "C:\Users\myname\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\ee_plugin\extlibs_windows\ee\data.py", line 340, in _execute_cloud_call
    raise _translate_cloud_exception(e) ee.ee_exception.EEException: Collection.loadTable: Collection asset 'users/myusername/master_shp/polygon' not found.


texts: Searching for changes of 2 layers based on geometry, Arcgis for desktop (ArcEditor

texts:
I´m trying to make a tool that would be able to compare 2 layers of any type (polygon, point, polyline) and find differences between them. Lets say i have a layer of roads (Polyline) from 2005 (layer X) and from 2015 (layer Y). This tool would check those layers and mark the roads missing in layer X (field Change = 1), roads missing in layer Y (field Change = -1) and the roads without change (field Change = 0). Note that matching objects between those layers by IDs might be uncorrect, I decided to make their own unique key.

So far I have a field named Geo_IDx and y (for both layers), that holds geometry token "SHAPE@JSON", created by arcpy.UpdateCursor:

    cursor2 = arcpy.da.UpdateCursor(Input2, [fieldnames2[0], "SHAPE@JSON", "Geo_IDx"])

Here comes the first issue: the string saved in "Geo_IDx" field is very long (currently set for 5000 signs, because I don't want to lose part of the geometry information).

Do you have any advice about how to make some better unique key for the comparison?

Another problem is that this comparison takes very long time for layers with big amount of objects stored in them.

Then I Merge the layers to have all those values in 1 layer, that holds the field "Change". After this I am comparing the "Geo_ID" fields from both layers, provided that if value from layer X matches with any value in layer Y, there was no change in that object. If the value is missing in either of those layers, something has changed, so by cursors I find in which layer the value is missing and update the Change value as mentioned above.

Following script shows the way I´m trying to update the "Change" field. geo1 and geo2 are lists of all values from Geo_IDx => geo1 and Geo_IDy => geo2.

    cursor3 = arcpy.da.UpdateCursor(Output,[f_names_merge[0], "Geo_IDx", "Geo_IDy", "Change"])
    for row in cursor3:
        if row[0]:
            if row[1] in geo2:
                row[3] = -1
            elif row[2] in geo1:
                row[3] = 1
            else:
                row[3] = 0
            cursor3.updateRow(row)

What I think the script is failing on is comparing the geometry strings, that are saved in Geo_ID fields. I dont have any idea why.

Questions:

Do you have any better idea of unique key for every object?
Is there any easier way to do this? (provided I´m using ArcEditor license).
Is it possible to compare the values of "SHAPE@JSON"?

I do not have access to FME.


texts: Applying single symbology definition to multiple TIFFs in ArcMap?

texts:
I have 143 separate TIFF files representing bathymetric data, all with different data ranges (high/low values). I want to create and apply a single symbology file having 6 levels corresponding to 6 discrete data ranges although individual TIFFs may only contain data with 1 or 2 of the 6 data ranges.

How do I do that in ArcMap 10.4?


texts: Dictate placement of InfoWindow - ArcGIS 3.27 API

texts:
Is there a way to dictate on which side of the point the InfoWindow shows up?

I have tried using map.infoWindow.show(location, InfoWindow.ANCHOR_UPPERRIGHT) thinking that this would make the window always show on the right side of the point, but it didn't seem to have an effect.

I'm trying to develop a way for the user to move the InfoWindow out of the way so that it doesn't block any other features that they may need to see, or in some cases block the feature that is highlighted. Even if I can't drag it out of the way, being able to move it 90 degrees would still help.


texts:
Is it possible to define a part of a geometry column as the Primary Key in a PostGIS database?

For example, when I want to store 3D points in a database, using the x/y coordinate as the primary key.

Thanks for helping!


texts: Why are two incompatible projections used together in OpenLayers?

texts:
While going through the examples in OpenLayers, I notice several examples with the following code:

map = new OpenLayers.Map('map_element',{
 ...
 units: 'm',
 projection: new OpenLayers.Projection('EPSG:900913'),
 displayProjection: new OpenLayers.Projection("EPSG:4326"),
});

My understanding is that the code:

projection: new OpenLayers.Projection('EPSG:900913')

displays a square-shaped map using a Spherical Mercator projection, and that the code:

displayProjection: new OpenLayers.Projection("EPSG:4326")

allows map controls such as the mouse to read from an ellipsoidal projection.

Notwithstanding the differences in units used (the former is in metres while the latter is in degrees), these two projections are totally incompatible since the points on these projections do not really match even if converted to the same unit.

So my question is: Why are these two projections being used together?


texts: QFIELD How to close the attribute table to be able to see the create icon again

texts:
I am trialing QField and when I click on a point it shows the table where I can view the attributes or edit them.

How do I close this window/tab so I can view the create feature pencil icon again?

I am sure it is obvious...


texts: Popups not showing on layer/feature/marker groups

texts:
I had two problems.  One was that, once I got lots of markers in one layer, and lots of polylines in another, someof my markers, starting from the north, started to vanish.
I noted in the doco that the zindex is derived initially from the latitude,so I bumped the offset up to +1000 for my markers, and -1000 for my polylines, and things look ok.
But I need the popups.
This suggests they actually dont work
https://github.com/Leaflet/Leaflet.markercluster/issues/180
So I've had to dump leaflets groups, and do my own layer management.
Further more, the popups seem to vanish when I use optional baselayers.
is this right ?


texts:
I'm trying to do a union between two Polygon Shapefiles, but when I do Vector > Geoprocessing Tools > Union I have the odd issue where most of it works fine, but some areas are returned in the union dataset as being a gap, but there is data there in both input files.

Here are some pictures to show you what I mean.

Input 1:

Input 2:

Union Output:

Union Output overlaid on the input layers:

To help identify the issue, they're both in the same CRS and are both Polygon format Shapefiles, and the geometry of the polygons in the input data where a gap is being created look fine under inspection.


texts:
I have created a popup that opens when the map is loaded. This is being used as a splash screen. I have updated the CSS with the following so it looks the part:

.custom-popup, .leaflet-popup-content-wrapper {
    background:#2c3e50;
    color:#fff;
    font-size:12px;
    line-height:10px;
    border-radius: 10px;
    }

.custom-popup, .leaflet-popup-content-wrapper a {
    color:rgba(255,255,255,0.1);
    }

.custom-popup, .leaflet-popup-tip-container {
    width:30px;
    height:15px;
    }

.custom-popup, .leaflet-popup-tip {
    background: transparent;
    border: none;
    box-shadow: none;
    }

JS:

var startpopup = L.popup({
    closeButton: true,
    autoClose: true,
    keepInView: true
    })
    .setLatLng(map.getBounds().getCenter())
    .setContent('
DISCLAIMER TEXT.
')
    .openOn(map);

However this is applied to all my popups that I have created for my different layers. I only want the above CSS applied to my splash screen popup.
How do I revert back to the original CSS for popups after the splash screen popup is removed?


texts: gdal2tiles.py generating dark tiles from .img file

texts:
Dark tiles are generated when I run gdal2tiles.py on a .IMG file(a scanned and georeferenced topographical map).
This question advised to apply a color relief map using gdaldem before generating tiles.
Does gdaldem work if the .IMG has no elevation data ?
I tried it using the color ramp below and the output file was all green.

0 110 220 110
900 240 250 160
1300 230 220 170
1900 220 220 220
2500 250 250 250

This is because the height for all raster cells are the same, 0.
Which color ramp should i use with gdaldem to prevent gdal2tiles.py from generating dark pngs ?
gdalinfo output is as below

Driver: HFA/Erdas Imagine Images (.img)
Files: ALOI.IMG
       ALOI.IMG.aux.xml
Size is 3830, 2552
Coordinate System is:
GEOGCS["GCS_WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["WGS_1984",6378137.0,298.257223563]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.017453292519943295]]
Origin = (32.999437318839924,3.000470692247897)
Pixel Size = (0.000392066994127,-0.000392066994127)
Metadata:
  DataType=Generic
Corner Coordinates:
Upper Left  (  32.9994373,   3.0004707) ( 32d59'57.97"E,  3d 0'1.69"N)
Lower Left  (  32.9994373,   1.9999157) ( 32d59'57.97"E,  1d59'59.70"N)
Upper Right (  34.5010539,   3.0004707) ( 34d30'3.79"E,  3d 0'1.69"N)
Lower Right (  34.5010539,   1.9999157) ( 34d30'3.79"E,  1d59'59.70"N)
Center      (  33.7502456,   2.5001932) ( 33d45'0.88"E,  2d30'0.70"N)
Band 1 Block=64x64 Type=UInt16, ColorInterp=Undefined
  Description = Layer_1
  Min=0.000 Max=255.000
  Minimum=0.000, Maximum=255.000, Mean=225.572, StdDev=31.338
  NoData Value=256
  Overviews: 1915x1276, 958x638, 479x319, 240x160
  Metadata:
    STATISTICS_COVARIANCES=982.046346716112,845.521289084631,710.510646703791
    STATISTICS_MINIMUM=0
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=225.57242379015
    STATISTICS_MEDIAN=0
    STATISTICS_MODE=0
    STATISTICS_STDDEV=31.337618714831
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_EXCLUDEDVALUES=
    LAYER_TYPE=athematic
  Image Structure Metadata:
    COMPRESSION=RLE
Band 2 Block=64x64 Type=UInt16, ColorInterp=Undefined
  Description = Layer_2
  Min=4.000 Max=255.000
  Minimum=4.000, Maximum=255.000, Mean=226.079, StdDev=29.739
  NoData Value=256
  Overviews: 1915x1276, 958x638, 479x319, 240x160
  Metadata:
    STATISTICS_COVARIANCES=845.521289084631,884.412989496947,765.062067015664
    STATISTICS_MINIMUM=4
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=226.07921779459
    STATISTICS_MEDIAN=0
    STATISTICS_MODE=0
    STATISTICS_STDDEV=29.739081853631
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_EXCLUDEDVALUES=
    LAYER_TYPE=athematic
  Image Structure Metadata:
    COMPRESSION=RLE
Band 3 Block=64x64 Type=UInt16, ColorInterp=Undefined
  Description = Layer_3
  Min=0.000 Max=255.000
  Minimum=0.000, Maximum=255.000, Mean=172.929, StdDev=27.135
  NoData Value=256
  Overviews: 1915x1276, 958x638, 479x319, 240x160
  Metadata:
    STATISTICS_COVARIANCES=710.510646703791,765.062067015664,736.335138826829
    STATISTICS_MINIMUM=0
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=172.92884721066
    STATISTICS_MEDIAN=0
    STATISTICS_MODE=0
    STATISTICS_STDDEV=27.135495920046
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_EXCLUDEDVALUES=
    LAYER_TYPE=athematic
  Image Structure Metadata:
    COMPRESSION=RLE

Is there another hack ?

Thanks,
Allan.


texts: How to create dynamic text of a layer name?

texts:
I'm trying to implement the code mentioned on this post:
Best way to have title which includes dynamic text of a layer name?
and I'm getting the error below. I am a novice and know nearly nothing about python.  I've been trying to troubleshoot on my own but this is over my head.

import arcpy, os

#Map document stuff
CurrentMXD = arcpy.mapping.MapDocument("CURRENT")
CurrentDF = CurrentMXD.activeDataFrame

#Output map directory
Output_dir = "C:\temp"

#Iterate through layer names
LayersToExport = ["Layer name 1","Layer name 2"]
for LayerName in LayersToExport:

#Get current layer, make visible
Layer = arcpy.mapping.ListLayers(CurrentMXD, LayerName, CurrentDF)[0]
Layer.visible = True

#Change title text
TextElement = arcpy.mapping.ListLayoutElements(CurrentMXD, "TEXT_ELEMENT","TitleText")[0]
TextElement.text = Layer.name

#Get output JPG filename, remove existing, export
Output_jpg = os.path.join(Output_dir,LayerName+".jpg")
if os.path.exists(Output_jpg): os.remove(Output_jpg)
arcpy.mapping.ExportToJPEG(CurrentMXD, Output_jpg, resolution=300, jpeg_quality=95)

#Make layer invisible again
Layer.visible = False

Runtime error
Traceback (most recent call last):
File "<string>", line 25, in <module>
File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\utils.py", line 182, in fn_
return fn(*args, **kw)
File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\mapping.py", line 1026, in ExportToJPEG
layout.exportToJPEG(*args)
File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\utils.py", line 182, in fn_
return fn(*args, **kw)
File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\arcobjects\mixins.py", line 568, in exportToJPEG
return self._arc_object.exportToJPEG(*args)
AttributeError: PageLayoutObject: Error in executing ExportToJPEG


texts:
Why is it that when I label my directories this way:

arcpy.env.overwriteOutput = True
arcpy.env.workspace = "C:\\Users\\Documents\\ArcGIS\\Testing"
arcpy.env.scratchWorkspace = "C:\\Users\Documents\\ArcGIS\\Testing"
scratch = "C:\\Users\Documents\\ArcGIS\\Testing"
space = "C:\\Users\\Documents\\ArcGIS" #General Directory
workspace = "C:\\Users\\Documents\\ArcGIS\\Textfiles" #Where the final product is saved
folder = "C:\\Users\\Documents\\ArcGIS\\AQS\\site_coordinates\\" #Location of the the site data
textfile_place = "C:\\Users\\Documents\\ArcGIS\\Textfiles\\pollution.txt" #Final Textfile
data = "C:\\Users\\Documents\\ArcGIS\\Textfiles\\data.txt" #Location of data

it reads all of my data, but when I try to clean it up and make it more transferable across computers like this:

arcpy.env.overwriteOutput = True
arcpy.env.workspace = "C:\\Users\\Documents\\ArcGIS\\Testing"
arcpy.env.scratchWorkspace = "C:\\Users\\Documents\\ArcGIS\\Testing"
space = "C:\\Users\\Documents\\ArcGIS"
scratch = space + "\\Testing"
workspace = space + "\\Textfiles"
textfile_place = workspace + "\\pollution.txt"
data = workspace + "data.txt"
folder = space + "\\AQS\\site_coordinates\\"

it doesn't work and says the files don't exist? I've been trying to look on the web but can't find a solution to the problem. Thanks!


texts: Using Cassini-Soldner projection in QGIS?

texts:
I have a shapefile in a Cassini Soldner reference system centered in Monte Mario (Italian cadastre map). I want to insert these data inside QGIS in order  to compare data with aerial images. When importing the file I have to select the reference system,  but I'm not sure  how to proceed.
I have noticed that the EPSG code for the Cassini Soldner representation is the 9806, but:

EPSG 9806 is not included in QGis, how I can import it?
Even if I successfully import it on QGis, how I can center it in Monte Mario Origin?

I am just a newbie about GIS, so please be understanding.


texts: Spatial adjustment of vector data in QGIS?

texts:
Is there any plugin or code reference available to spatially adjust vector data based on points selection like in the Spatial Adjustment tool of ArcGIS Desktop?

I am trying to develop spatial adjustment tool plugin in QGIS.


texts: How to change the active layer in QGIS Console?

texts:
I want my script to perform actions on a particular layer.  Online advice I've found so far suggests using the iface method "setActiveLayer()".  Makes sense based on the name alone.

Scenario:  In the Layers window I have a layer named "Province" highlighted (active).  Then I run the following code:

>>> vl = QgsMapLayerRegistry.instance().mapLayer( 'na_roads_Prov' )
>>> iface.setActiveLayer(vl)

Console Output = False

As far as I understand executing this code in the console should move the highlighted "active layer" from "Province" to 'na_roads_Prov'.  Yes?  Instead the "Province" layer remains highlighted as active, and the console output says "False".  Any ideas?  Am I not referring to the layer correctly?


texts: Script fails when executed from toolbox but not in command window (MakeFeatureLayer)

texts:
I am able to run the below code in the Python command window in ArcMap 10.2.1 without issue. However, when I run this as a toolbox script, it fails at line 15:

  newLyr = arcpy.mapping.ListLayers(mxd, "New Feature Layer", df)[0]

  IndexError: list index out of range

It seems to me that the feature layer is not being generated in line 14, causing the error in the following line. But I'm not sure. Does anyone have any idea why the processing would fail when executed as a script but not when run in the command window?

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
arcpy.Delete_management("New Feature Layer") #clear any earlier instance from memory

attrLyr = arcpy.mapping.ListLayers(mxd, "SELECT FEATURES BY ATTR LAYER", df)[0]
locLyr = arcpy.mapping.ListLayers(mxd, "SELECT FEATURES BY LOCATION LAYER", df)[0]
srcLyr = arcpy.mapping.ListLayers(mxd, "SRC SYMBOLOGY LAYER", df)[0]

expression = "SELECTION EXPRESSION"

arcpy.SelectLayerByAttribute_management(attrLyr,"NEW_SELECTION",expression) #select features from attrLyr according to expression
arcpy.SelectLayerByLocation_management(locLyr, "INTERSECT", attrLyr) #select features from locLyr that intersect with attrLyr
arcpy.MakeFeatureLayer_management(locLyr, "New Feature Layer") #create feature layer from selected features in locLyr
newLyr = arcpy.mapping.ListLayers(mxd, "New Feature Layer", df)[0] #assign new layer to variable "newLyr"
arcpy.mapping.UpdateLayer(df, newLyr, srcLyr, True) #apply srcLyr symbology to newLyr


texts: Contour lines with imported shapefile using Blender GIS - wrong Z values

texts:
I've been making experiments with Blender GIS and QGIS.

I'm having problems when importing Shapefiles with contour lines, with a field for elevation (ELEVATION FROM FIELD selected when importing).

I'm describing them below (Blender 2.75 and QGIS 2.10).

When there are contour lines that have no value for elevation (in QGIS for example, Z is NULL) mixed with contour lines that have values, Blender GIS cannot load the SHP. This is no real problem actually, since it's easy to fix in QGIS.
Now the real problem: when I import a shapefile into Blender with all the contour lines with Z values, it worked with a shapefile with only 3 contour lines (3 features in QGIS), with Z values 20, 40 and 60.

But when I imported a shapefile with 49 contour lines (49 features in QGIS) with Z values from 0 to 800 (20m interval), the import result was very messy. It appeared that the contour lines assumed random values, but still inside the original values from the shapefile.

Now the fun part: when importing, instead of selecting ELEVATION FROM FIELD, if I select EXTRUSION FROM FIELD, everything is right.

Then I can just delete the 0 value vertices and I end up with just the contour lines in the right Z position (so I can apply Delaunay triangulation to have a surface).

So, conclusion is, I can reach the final result I want, but it needs a bit more working.

Anyone have an idea of what is going on with ELEVATION FROM FIELD selected?


texts:
I managed to load osm-data with osm2po, pgrouting, postgis and postgresql into qgis.

Then I did a drivingdistance analysis based on cost to finally get the qgis contours plugin up and running.

But when I run it, it doesnt really output "the truth".

I want to get also the cost max-points (in my case 1), but somehow the contours are ending somewhere at 0.7?

What am I doing wrong? please see image attached.


texts: Which world map projection allows comparing countries by shape and size?

texts: The Mercator projection increases scale along with the latitude. Greenland seems larger than Australia, when in fact it is much smaller. The equal-area projections however stretch shapes away from the Equator. I suppose that there exist no projection that preserves size and shape at the same time. But for the purpose of comparing countries it is in fact not needed to preserve any properties exactly. The distortions should be only roughly the same in any point of the map. What is the map projection that accomplishes it most closely?


texts: Eliminating single raster pixels of flow distance raster in ArcGIS Desktop?

texts:
I want to clean up a calculated flow distance raster and need some ideas on how to do it.

As you can see in the picture, there are some pixels along streets and rivers (gray) that show connection (Red circle). These are the ones I want to get rid off.

Single pixels further away from the streets or rivers should not get deleted.
Pixels close to streets/rivers within the bigger cluster areas should not be deleted.

Therefore, I am looking for a way to select pixels with only 2 neighbor pixels along the streets and delete them.

I'm thinking about a workflow like this:

Select all Pixels close to the streets and rivers
Calculate how many pixels are close to each pixel (White are NoData)
Delete if pixel only touches 2 other pixels


texts: Which version of Sen2Cor is used by the ESA?

texts: I'm trying to find out which version of the Sen2Cor processor was used by the ESA to process a Sentinel L2A product for a specific date. The latest version is Sen2Cor v2.8 which was published this year. However, for data of the last two years the ESA had to use a different version. I can't find any information about the used Sen2Cor in the metadata when I am download L2A products via the Copernicus Open Access Hub. This information would be necessary to decide whether to process a Sentinel L1C product with the new Sen2Cor processor or not.


texts: Extracting data from a .mat file with numpy and a mask

texts:
I have the output of a xBeach model run in .mat format. There are several layers that represent the sea status for each time-step. In each layer there are the coordinates and the water level of the relative point.
I've succeeded in opening it through Python with scipy and, for example, requesting the maximum of each layer.
    import scipy.io
    import os
    import numpy

    FILE = 'level.mat'

    layers = scipy.io.loadmat(FILE)['data'][0][0][0]

    print(*map(numpy.nanmax, scipy.io.loadmat('FILE')['data'][0][0][0]))

Now I'd need to select the maximum of each layer, but in a specific area. To determine this area I've a polygon as geojson or a shapefile that I could use, but I still have troubles in finding how to do to combine this with the layers file as numpy array and extract the maximum in that area.
Does anyone have an idea?


texts: Conversion between WGS84 height and local height datum

texts:
For x/y or lat/lon or heights between different ellipsoid, we can do conversion by using proj.4. But I found that if I want to convert between WGS84 height and local height, I found that it is not so simple.

For example, I try to convert height of EPSG:5738 which is Hong Kong Principal Datum from WGS84 height. As it proj.4 file is blank, I cannot do it.

From the official description of HKPD (http://www.geodetic.gov.hk/smo/gsi/Data/pdf/explanatorynotes.pdf), it just shows the rough relationship.

Does it mean I have no luck with doing it in proj.4 or other tools?


texts:
Is it possible to read a shapefile with ogr as a numpy array, in order to use it as a mask to filter out a geotiff image?

What I would like to obtain is a numpy array from the layer below:

from osgeo import ogr

map = ogr.Open('file.shp')
layer = map.GetLayer()

When I rasterize the shapefile using gdal_rasterize (into tif), I cannot find the correct resolution for the obtained raster that matches the resolution of my geotiff file.

Note: The shapefile I'm trying to read contains around 6K features with 5 fields, among them is the area_ha.

gdal file.tif (a subset of a Sentinel-1 SAR image):

Driver: GTiff/GeoTIFF
Files: preimage.tif
Size is 9146, 11277
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (97.254882569202280,17.422529922059596)
Pixel Size = (0.000089831528412,-0.000089831528412)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_IMAGEDESCRIPTION=preimage
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  97.2548826,  17.4225299) ( 97d15'17.58"E, 17d25'21.11"N)
Lower Left  (  97.2548826,  16.4094998) ( 97d15'17.58"E, 16d24'34.20"N)
Upper Right (  98.0764817,  17.4225299) ( 98d 4'35.33"E, 17d25'21.11"N)
Lower Right (  98.0764817,  16.4094998) ( 98d 4'35.33"E, 16d24'34.20"N)
Center      (  97.6656821,  16.9160148) ( 97d39'56.46"E, 16d54'57.65"N)
Band 1 Block=9146x11277 Type=Float32, ColorInterp=Gray

ogrinfo map.shp:

INFO: Open of `ST20150315_Preflood.shp'
      using driver `ESRI Shapefile' successful.
1: ST20150315_Preflood (Polygon)


texts: Problem in importing shapefile to Postgis with postgis Shapefile manager

texts:
i have a problem with importing my shapefile to postgreSQL, it's a poroblem with the encoding, here's the screenshot:

i tried to change encodingin "options" buti get this message:

"Failed in pgui_exec(): ERREUR:  séquence d'octets invalide pour l'encodage « UTF8 » : 0xe3 0xa9 0x72

Shapefile import failed.
"

Do i have to change the encoding in arcgis?and how?

thnks in advance


texts:
Ok so I've got a noodle scratcher and want some input before I open my mouth to my higher-ups.  I recieved a new boundary description from a surveyor that uses metes and bounds in the form of bearings N XX°XX'XX" W going XX feet and so on.  The coordinate system they used was a custom Traverse Mercator projection.  However, my state, North Dakota, uses State Plane S Nad83 Int. Feet (a Lambart Conformal Conic projection) as its standard coordinate system; so I believe that all their records would be stored in that system.

My question is what, if any, difference will there be in the bearings if it were transformed into the State Plane system?  Both systems use the same units, datum, and spheroid.

I understand their reasoning for using the custom coordinate system as they shifted the false easting and northing to reduce distortion in the location they preformed the survey and using a perfect square coordinate system without converging lines of longitude make sense at this scale.  However, because all the records I see for the state never specify a coordinate system I think it is implied that the state's standard coordinate system is used.  I don't want this data on file without the custom coordinate system attached and have someone in the future go out and use the description using the wrong coordinate system.  Or is the discrepancy between the two, if any, so miniscule that I'm being needlessly nit-picky?


texts: Tips and advice for editing complex Models / Tools in ModelBuilder ArcGIS - Advanced issues

texts:
I am building a series of Modles / Tools in ArcGIS ModelBuilder and so I am frequently editing different Models / Tools within Toolboxes. The Tools generally have at least two submodels and contain upwards of 30 separate Tools / processes.

I am using standard Modelbuilder Tools and Toolboxes, not Python script versions.

The issues I am requesting answers on is the behavior of models / Tools that begin to appear unstable and start throwing up Errors that are inconsistent between model runs or that are frequently solved by closing down ArcCatalog and then re running the model / Tool. It almost appears that some models after lots of editing or cutting and pasting sections of  models between different models become "unstable" and throw up errors.

My question is this - do people have advice or help on how to maintain the stability of models / Tools in ModelBuilder when they are being edited for long period of time, and becoming very complex with lots of nested submodels?

Is there a best praqctice for tracking different versions of Model, e.g. renaming etc.

I am using ArcGIS Desktop Advanced 10.1

I have not posted a model / tool example because the issues i am thinking of have occurred sporadically across several different models and I have checked through all possible normal errors resulting from error codes etc.

The types of errors are:
- model "unable to save" or "failed to save"
- model becomes unstable and crashes, when sub models are opened
- field calculator using VBA code suddenly stops working, and only works when changed to Python code version even though previous runs of the model have worked fine using the VBA field calculator
- workspace or data source is read only

My current solutions to these issues have been:

use in_memory for calculation where possible and keep deleting files that build up in  memory after they have been used
setting scratch and inputs, outputs work spaces within the model as parameters
emptying the scratch workspace when lots of files build up from multiple model runs
copy and paste the model into a new empty model and save

Any thoughts, help or tips from advanced ModelBuilder users would be greatly appreciated.

As this question was closed for being too vague I have opened some related questions here to hopefully share ModelBuider expertise and best practice:

Do you use a proven workflow to edit and document ArcGIS ModelBuilder models?

and

Troubleshooting why ArcGIS Model works on one PC but not another?

and

Avoiding or dealing with "unable to save" errors in ArcGIS ModelBuilder?


texts: Generating STL (3D) file from raster ASTER GDEM data?

texts:
I load an ASTER data to Global Mapper workspace and try to export to STL file, receive the message that I need a 3D vector file.

I tried all formats in GM, grid xyz, dem, contour in vector format and so on.

What kind of format file I need to convert ASTER data to export to STL file?

or even, PLY, OBJ,SKP, 3D formats?


texts: Adding system tools to custom toolbox in ArcGIS Pro?

texts:
I am still getting over the quirks of ArcGIS Pro, right now I want to create my own toolbox that contains some tools useful for my current project: hillshade, contour and kernel density.

But it seems to me that there is no way to add system tools to a custom toolbox.

Can this be done?


texts: Speed difference MS SQL Server 2019 to PostgreSQL/PostGIS 12 for spatial data

texts:
I imported a shapefile with 3000 polygons in UTM 25832 into both databases MS SQL Server 2019 and PostgreSQL/PostGIS. Then I calculated a buffer of 1000 meters around these geometries and transferred the result to a new table. (Again in both databases with their integrated STBuffer() functions)
Result:
PostGIS: 9 seconds
MS SQL 2019: 1 minute
MS SQL 2017: 1:09 minutes
All databases are running on the same PC and both MSSQL Server and PostgreSQL/PostGIS have been freshly set up. No optimizations were deliberately carried out and no spatial indices were used.
Is it really possible that PostgreSQL/PostGIS is so much faster than MS SQL Server 2019 as far as geodata is concerned?


texts: Messages in JSON output for a Geoprocessing Service from a Model

texts:
I have a Model created from ModelBuilder. It was published as a Geoprocessing Service on ArcGIS Server 10.1

When it was triggered from the REST service, the JSON output printed out the "results" and "messages". The "results" can be configured as derived parameters in ModelBuilder. However, I  couldn't find anywhere that "messages" can be configured.

Is it possible?


texts:
I am looking for a way to create a 10 second graticule/grid. I have north west corner coordinates where the grid should start from and spread out based on 10 second intervals. I have looked at two arcmap tools but they didn't work for me. The "Grid Index Features" tool worked well except for the fact that it places my starting coordinates at the south west corner and I don't see an option to have it start in the north west. The "Create Fishnet" tool parameters are shown below but they result in a fishnet which starts in the north east corner instead of the north west corner. I used Nad83 and converted 10 seconds to decimal degrees and used that as the cell size width and height.

The green point is where I need the fishnet to start from and draw to the east but it by default, starts at the north east corner and draws to the west.

Question 1. Is there a way to have the fishnet start drawing from the north west corner spreading east?
Question 2. How can I verify that the intervals are indeed 10 seconds apart?


texts: Draw lines between points from two layers in QGIS

texts:
In QGIS 3.10.4 on Mac, I have two multipoint shapefiles: one is "Start" points, and the other is "To" points. Each layer has "x" and "y" coordinate columns in the attribute table. "Start" points only have "location_name" (string) as IDs, and "To" points have ID numbers as IDs.
I also merged the two layers, so I have a layer that has both sets of points if needed.
My goal is to draw lines between each "start" point and each "to" point. Many plugins in other answers cannot be found now.
I'm wondering how can I do so?
The points look like the picture below, red dots are Start points, and blue ones are Tos.


texts: Overpass API: Query points not inside polygons

texts:
I need to query the Overpass API for some data that comes down to the question:

  Give the set of nodes that is not contained (in spatial terms) in a set of areas (given by closed ways).

I looked into Overpass documentation and on the wiki (https://wiki.openstreetmap.org/wiki/Overpass_API/Overpass_QL#Map_way/relation_to_area_(map_to_area) and https://wiki.openstreetmap.org/wiki/Overpass_API/Overpass_API_by_Example), but I fail to see how to use area, and other query methods can make my query.

The areas (residential areas): https://overpass-turbo.eu/s/SgN

The points (places): https://overpass-turbo.eu/s/SgQ

My best try that doesnt work (https://overpass-turbo.eu/s/SgO):

area["landuse"="residential"]({{bbox}})->.searchArea;

node["place"](area.searchArea)->.a;
node["place"]({{bbox}})->.b;

(.b; - .a;);

The entire query can be scoped to bbox.


texts: Customizing User roles in Portal for ArcGIS

texts:
Portal for ArcGIS allows creating custom roles with different privileges defined. Fine. But if I need even more detailed customization?

For example if I need to create a role which can create a map but not an app? Is there a way to separate the item "Content" to Maps and Apps and define different privileges for each?


texts: How to know what GRASS version are your python scripts using?

texts:
When you are running python scripts in QGIS via the QGIS own's python console or from a QGIS plugin (when using algorithms from the processing framework).

How to know or even change the version of GRASS that is being used?

This question came up while discussing this other question.


texts: Limit Nominatim results tp specific country

texts:
I want to limit the search in Nominatim to a specific country and was trying this code:

 mapPanel = new GeoExt.MapPanel({
    border: true,
    region: "center",
    map: map,
    zoom: 7,
    layers: [
            locationLayer = new OpenLayers.Layer.Vector("Addresses", {displayInLayerSwitcher: false,
            styleMap: new OpenLayers.StyleMap(new OpenLayers.Style({}, {rules: rules}))
                })
    ],
    tbar: [topbarItems],
     bbar: [{
        xtype: "gx_geocodercombo",
        layer: locationLayer,
        hideTrigger: true,
        emptyText: "Search address",
        zoom: 10,
        minChars: 4,
        queryDelay: 100,
        countrycodes: "de",
        srs: "EPSG:4326",
        width: 300
    }
     ]

});

The wiki has a parameter listed that goes by the name: countrycodes and later on the same page there is another one called: country_code. It doesn't matter what parameter I use in the code above, when typing in the combobox, the results are coming in from all over the world.

Can somebody here tell me on how to limit the search to a specific country?


texts: Rectangle (bbox) around non masked values in GEE

texts:
I am looking for a solution to create a bounding box (rectangle) around non-masked pixels of an image in Google Earth Engine.
In this simple example, the first image has only some non-masked pixels within the larger region of interest (roi). I would like ot automatically create a rectangle that includes all of the existing values within the roi.
// cloud maks
function cloudMask(image){
  var QA60 = image.select('QA60');
  return image.updateMask(QA60.lt(1));
}

var S2dat = ee.ImageCollection("COPERNICUS/S2")
        .filterDate('2018-06-01', '2018-06-15')
        .filterBounds(roi)
        .map(cloudMask)
        .first();

Map.setCenter(126.42, 72.80, 6);
Map.addLayer(S2dat.select("B2").clip(roi));

https://code.earthengine.google.com/c8e8a6b86f8fb8023151f4c145a36901


texts: Controling affine transform of raster in ArcGIS Pro

texts:
I have a png file which I created from a field equation. I need to plot that on my map. I am plotting it on a map with Coordinate System being GCS WGS 1984 (EPSG: 4326).

Since it is not georeferenced, the image comes in with the top left corner at N 0 E 0, which I anticipated. It also comes in with pixel size at 1 degree. I need to affine transform it by specific numbers, and not by clicking control points.

I have looked around online for a way to do this, but I can't find anything. Can someone point me in the right direction?

If necessary, I can do this in ArcGIS Desktop 10 then move it to ArcGIS Pro.


texts: Where to download free transport and urban data for the capital cities in South East Asia (ASEAN)?

texts:
I'm currently looking for the transport (mainly public transport - rout lines and stop locations) and urban (public amenities, population density etc) data on the capital cities in South East Asia (ASEAN Region).
Any idea where I might find them?


texts: OpenLayers (allOverlays): Vector layer not rendered without raster layer beneath it

texts:
This is another newbie question. I'm trying to create a map using OpenLayers that comprises only of vector layers from a WFS server. However, I'm thwarted because my vector layer won't render unless it has a raster layer beneath it. I surmise this is because OpenLayers is deriving some setting(s) (resolution? minExtent?) from the raster layer that it's unable to do from the vector one, but I'm just guessing and I don't know what those settings are.

The code below works beautifully, displaying both layers. If I remove the line to add the raster layer to the map, or add it after I've added the vector layer, the map displays nothing (calling oMap.zoomToMaxExtent() still shows nothing).

What am I missing?

function initMap() {
    var oRasterLayer,
        oVectorLayer,
        oBounds,
        oSpatialProj,
        oSphericalMercatorProj;

    // Use MapBox dark theme for UI controls
    OpenLayers.ImgPath = "http://js.mapbox.com/theme/dark/";

    oSpatialProj           = new OpenLayers.Projection("EPSG:4326");
    oSphericalMercatorProj = new OpenLayers.Projection("EPSG:3857");

    oBounds = new OpenLayers.Bounds(-11, 49.6, 4, 62.6).transform(
        oSpatialProj,
        oSphericalMercatorProj
    );

    oMap = new OpenLayers.Map(
        sMapDivId,
        {
            numZoomLevels    : 10,
            maxExtent        : oBounds,
            restrictedExtent : oBounds,
            allOverlays      : true
        }
    );
    oMap.addControl(new OpenLayers.Control.LayerSwitcher());

    oRasterLayer = new OpenLayers.Layer.XYZ(
        'Land',
        sTilesUrl.replace(/\${layer}/, 'land'),
        {
            projection        : oSphericalMercatorProj,
            wrapDateLine      : true,
            sphericalMercator : true
        }
    );

    /* Add a postcode area selection layer. This is a vector layer, working
    directly with data from PostGIS. */
    oVectorLayer = new OpenLayers.Layer.Vector(
        'Selection',
        {
            projection        : oSphericalMercatorProj,
            wrapDateLine      : true,
            sphericalMercator : true,
            strategies        : [new OpenLayers.Strategy.BBOX()],
            protocol          : new OpenLayers.Protocol.WFS({
                url          : '/wfs',
                featureType  : 'areas',
                featureNS    : "http://gis.mediatel.co.uk",
                geometryName : 'wkb_geometry',
                srsName      : 'EPSG:3857'
            })
        }
    );

    // SWAP THESE TWO LINES, OR COMMENT-OUT THE RASTER ONE, TO BREAK THE MAP:
    oMap.addLayer(oRasterLayer);
    oMap.addLayer(oVectorLayer);

    oSelectControl = new OpenLayers.Control.SelectFeature(
        oVectorLayer,
        {
            clickout    : true,
            toggle      : true,
            multiple    : true,
            hover       : false,
            toggleKey   : 'ctrlKey',
            multipleKey : 'shiftKey',
            box         : false
        }
    );
    oMap.addControl(oSelectControl);

    oMap.setCenter( // Near London
        new OpenLayers.LonLat(0, 51.5).transform(
            oSpatialProj,
            oSphericalMercatorProj
        ),
        9
    );
}


texts: what is the solution for this funky sorting Esri does in attribute tables with decimal places values... screen shot: http://i.imgur.com/0fqvF.png attribute table funky sort


texts: Getting different resolution between horizontal/vertical directions, from lon/lat to openlayers map?

texts:
I draw several points of Lon/Lat(EPSG:4326) on the OpenStreetMap v 4.6.5 (EPSG:3857).

It seems like the distance between the points are closer in the Latitudinal direction then in the Longitudinal direction.

The resolution/distance is approx 5 cm between the horizontal points, and approx 10 cm between the longitudinal points (Measured with the OSM-measure plugin).

Near equator this does not seem to be a problem, as the distance between all points seems uniform.
I thought plotting lon/lat points on OSM, should be straight forward :)

What is causing this problem and how do i solve this further north of equator?


texts: Creating Python script to reconcile/post geodatabase versions?

texts:
I’m trying to set up a python script so that I can schedule our rec/post to run on a server.  I’ve created a model with the Reconcile Version tool, with the option to Post checked on. I then export this to a Python script.  I run the script, and it finishes without any errors, but I do not see that my edits are posted to the default version. Is it not possible to script this tool, or is there a switch that I’m missing?

I’m using ArcSDE 9.3.1.


texts: Hi I'm working with OpenLayers version 3 and using method map.forEachLayerAtPixel(). The problem is I usually work on big screen and now I have switched to my laptop with smaller one and suddenly function is not working at certain zoom levels. I have to zoom at max and only then it works. Does someone encountered this problem maybe ? The same method when working with OpenLayers5 works properly..


texts: Seeking Brazil wide reference system to compute road length?

texts:
I've a Brazil road layer that include start km and end km for each road, some are empties.
Have checked and some of the values are wrong so I'll calculate road length by myself. I use ArcGis, but can use any other (Qgis, postgresql-postgis).

Original data is in SIRGAS_2000
GCS_SIRGAS_2000
WKID: 4674 Authority: EPSG

Also have available:
SIRGAS_2000_Brazil_Mercator
WKID: 5641 Authority: EPSG

What's Best Brazil wide reference system to compute road length?


texts: Why are geolocations on Twitter "rounded"?

texts:
I'm plotting some geolocations (lon, lat) of tweets collected using the Stream API with a 7mi. radius. The picture below shows two maps of the same data, but using different alpha/transparency for each observation.

Looking at the map on the left, I noticed that the geolocations follow mostly a sort of imaginary "grid" with specific nodes.
Anyone could help me to understand why this pattern surfaces?

Of course, I do have roundings in the coordinates, eg:

"loc: 42.7388,13.1798"         "loc: 42.6252,13.2948"         "loc: 42.6008,13.293"          "loc: 42.73,13.2028"
"loc: 42.66918468,13.27893702"

but my question is why these roundings occur since I assumed that the Twitter API would guarantee the same level of precision?


texts: QGIS Project saved in PostgreSQL not loading on other computers

texts: I have created QGIS projects that are reading from PostgreSQL(data source). I then saved the individual projects to PostgreSQL in the public schema under their respective databases. When I open or load the projects on my machine from PostgreSQL, they open fine without any issue. However, when I open the projects from PostgreSQL on other computers I cannot see them. I cannot see the icon when logged on other computers whether as a superuser (I can only see the rest of the other layers). Also when I log in as a defined user with read-only/view rights I cannot see the projects.


texts: Exporting region of orthomosaic using QGIS to TIFF lacks TFW file

texts:
I have a very large orthomosaic that's come via Pix4d, which I'm trying to move into AutoCAD. Because its so large I understand that Pix4d use their own proprietary .TIF format which AutoCAD can't read. QGIS is happy with it so I've been dragging it over to QGIS and then exporting regions.
Recently when I've tried to export a region it doesn't save a separate .TFW file. Just the TIFF. How do I force this?


texts: Creating Centrelines from Road Polygons/Casings using ArcGIS Desktop?

texts:
I want to create center lines from road polygons.

Are there any automated tools in ArcGIS for Desktop to do this?

I am including a Screenshot of my vector file below.


texts: Geoserver: Increase the precision of the json return of WPS call

texts:
I'm calling the geo:centroid WPS process from Geoserver. If I set the output type to "application/json" the returned values has only 4 digits after the dot (e.g xx.xxxx), but if I set the output type to "text/xml" the returned values has more digits after the dot (e.g xx.xxxxxxxxxxxxxxx).
My question is how to make the "application/json" returned values have more precisions (more digits after the dot).

Here is an example:

1- application/json

<?xml version="1.0" encoding="UTF-8"?>
<wps:Execute version="1.0.0" service="WPS" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.opengis.net/wps/1.0.0" xmlns:wfs="http://www.opengis.net/wfs" xmlns:wps="http://www.opengis.net/wps/1.0.0" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:wcs="http://www.opengis.net/wcs/1.1.1" xmlns:xlink="http://www.w3.org/1999/xlink" xsi:schemaLocation="http://www.opengis.net/wps/1.0.0 http://schemas.opengis.net/wps/1.0.0/wpsAll.xsd">
  <ows:Identifier>geo:centroid</ows:Identifier>
  <wps:DataInputs>
    <wps:Input>
      <ows:Identifier>geom</ows:Identifier>
      <wps:Data>
        <wps:ComplexData mimeType="application/json">
            <![CDATA[
        {
  "type": "LineString",
  "coordinates": [
    [
      [
        0,
        0
      ],
      [
        10.25124,
        0.257741
      ],
      [
        10.578814,
        10.255447
      ]
    ]
  ]
}
]]>
        </wps:ComplexData>
      </wps:Data>
    </wps:Input>
  </wps:DataInputs>
  <wps:ResponseForm>
    <wps:RawDataOutput mimeType="application/json">
      <ows:Identifier>result</ows:Identifier>
    </wps:RawDataOutput>
  </wps:ResponseForm>
</wps:Execute>

it returns

{
  "type": "Point",
  "coordinates": [
    7.7375,
    2.6609
  ]
}

2- text/xml

<?xml version="1.0" encoding="UTF-8"?>
<wps:Execute version="1.0.0" service="WPS" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.opengis.net/wps/1.0.0" xmlns:wfs="http://www.opengis.net/wfs" xmlns:wps="http://www.opengis.net/wps/1.0.0" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:wcs="http://www.opengis.net/wcs/1.1.1" xmlns:xlink="http://www.w3.org/1999/xlink" xsi:schemaLocation="http://www.opengis.net/wps/1.0.0 http://schemas.opengis.net/wps/1.0.0/wpsAll.xsd">
  <ows:Identifier>geo:centroid</ows:Identifier>
  <wps:DataInputs>
    <wps:Input>
      <ows:Identifier>geom</ows:Identifier>
      <wps:Data>
        <wps:ComplexData mimeType="application/json">
            <![CDATA[
        {
  "type": "LineString",
  "coordinates": [
    [
      [
        0,
        0
      ],
      [
        10.25124,
        0.257741
      ],
      [
        10.578814,
        10.255447
      ]
    ]
  ]
}
]]>
        </wps:ComplexData>
      </wps:Data>
    </wps:Input>
  </wps:DataInputs>
  <wps:ResponseForm>
    <wps:RawDataOutput mimeType="text/xml">
      <ows:Identifier>result</ows:Identifier>
    </wps:RawDataOutput>
  </wps:ResponseForm>
</wps:Execute>

it returns

<?xml version="1.0" encoding="UTF-8"?>
<gml:Point xmlns:gml="http://www.opengis.net/gml" srsDimension="2" xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:sch="http://www.ascc.net/xml/schematron" xmlns:xlink="http://www.w3.org/1999/xlink">
    <gml:pos>
        7.737501114169446 2.6609131596654136
    </gml:pos>
</gml:Point>


texts: Converting SpatialLite.Core.API.IGeometry to GeoAPI.Geometries.IGeometry

texts:
I have an in-house .NET library which processes geometry data in GeoAPI form (https://github.com/NetTopologySuite/GeoAPI). Within this library I have a method which I need to call which is expecting a GeoAPI.Geometries.IGeometry parameter.

I'm processing OSM PBF data using the PbfReader in SpatialLite.Osm (http://spatial.litesolutions.net/index.html). This gives me objects which implement SpatialLite.Core.API.IGeometry.

My question is how can acheive a conversion from SpatialLite.Core.API.IGeometry to GeoAPI.Geometries.IGeometry? I've made a few searches and found nothing to help me. I can go about creating a conversion procedure myself but it strikes me that there might be a better way.

Any suggestions?


texts: Recursively getting all lyr file properties from group layer on disk using ArcPy?

texts:
I'm writing a python script tool that will recursively search through a folder and all subfolders and find all .lyr files.  Ultimately, I'm writing a CSV file which contains the path of the .lyr file, the name, the source, and the definition query.

So far, so good.  I'm able to do this with no troubles.

The tricky bit is to get the source of all sub-layers in a group layer.  Because group layers can be nested and sub-nested, etc, I decided to write a recursive function to get the properties of the layers.  For the sake of clarity in this question, I've modified the code to just print instead of returning all the properties of the actual function.

Code for the recursive function:

def searchLayer(layerObj):
    if not layerObj.isGroupLayer:
        print layerObj.longName
    else:
        for sublayer in layerObj:
            searchLayer(sublayer)

Again, this works for the most part.  However, it's not printing the longName property as I would expect.  Suppose I have a .lyr file called "Administrative Boundaries".  And it looks like this:

Administrative Boundaries
    - Province Boundaries
    - Local Government Boundaries
    - National Boundaries

What I would expect to see as an output is the following:

Administrative Boundaries\Province Boundaries
Administrative Boundaries\Local Government Boundaries
Administrative Boundaries\National Boundaries

However, what I see instead is:

Province Boundaries
Local Government Boundaries
National Boundaries

As you can see, it's not printing the full longName path as it should.  What am I doing wrong here?

I'm using ArcGIS Desktop 10.3.1 - Basic License.

If I run the following code in an MXD, I get the expected results.

mxd = arcpy.mapping.MapDocument('current')
layerList = arcpy.mapping.ListLayers(mxd)
for layer in layerList:
    print layer.longName


texts: How to disable layer uploading in GeoNode?

texts: I'm using GeoNode on my Ubuntu 16.04. I want to disable layer publishing for all users and enable that just for a group of users. How can I do that?


texts: How to handle licenses when reinstalling ArcGIS 10.1 after computercrash?

texts:
In our student organisation we are having some issues with ArcGIS 10.1. When a computer crashes with ArcGIS 10.1 and with a valid student license for one year, is it possible to reinstall it after the crash with the same license number? A lot of our students had some real troubles with it because it says it's already installed on a machine, and the only thing our university gives as an answer is: "just come and swap your license with a new one".

Even in the ArcGIS Support Center I can't find a good solution for this.
Hopefully someone here got a solution for this or know a way-around for the problem, because a lot of valid licenses were already thrown away.


texts:
In PyQGIS, I kwnow how to draw polygons starting from a list of points which represent the nodes of the shape. An example of list could be this:

# First polygon
list_1 = [QgsPoint(0,0),QgsPoint(0,50),QgsPoint(50,50),QgsPoint(50,0)]

# Second polygon
list_2 = [QgsPoint(10,10),QgsPoint(10,30),QgsPoint(30,30),QgsPoint(30,20)]

which lead to the representation of these polygons:

If I manually do some operations, I can create a unique polygon like this one:

My question is: Is there a specific rule or a particular way to handle points (in or outside a list) for obtaining the result above? I only know that, for simple polygons, the points must to be ordered clockwise.

My goal is similar to that one in this question, but I think that the context is different here.


texts: Hide 'geometry' object on WFS GetFeature request for multiple TypeNames

texts:
I'm trying to fetch features for a WFS - but I'm not interested in the Geometry object. I would like to exclude it from the results...

Is there any way to do so?

I managed to define the TYPENAME parameter with the list of properties I needed for 1 specific Type. Following does not return the Geometry object. Perfect for this case...

SERVICE=WFS
&REQUEST=GetFeature
&VERSION=1.1.0
&MAXFEATURES=9999
&CRS=EPSG:31370
&BBOX=142562.4355,185629.719,142722.4355,185789.719
&OUTPUTFORMAT=application/json
&TYPENAME=Type1
&PROPERTYNAME=Property1,Property2

Problem

I want to call multiple types in a single GetFeature request.

I specify the list of TYPENAME...
I also define the list of PROPERTYNAME...

But the list of PROPERTYNAME contains values which don't occure in all TYPENAME items.

SERVICE=WFS
&REQUEST=GetFeature
&VERSION=1.1.0
&MAXFEATURES=9999
&CRS=EPSG:31370
&BBOX=142562.4355,185629.719,142722.4355,185789.719
&OUTPUTFORMAT=application/json
&TYPENAME=Type1,Type2,Type3,...
&PROPERTYNAME=Property1,Property2,Property3,...

Response as Exception:

Requested property:  Property3 is not available for Type1. The possible propertyName values are: ...

Following does not seem to change anything

&ows/wfs_geomtype=None

Is there any way to exclude the GEOMTRY object with above WFS GetFeature request?


texts:
I know this has been asked before, and there have been a lot of responses to this;

https://community.esri.com/t5/arcgis-pro-questions/why-does-arcgis-pro-have-to-be-so-slow/td-p/599220
Redrawing slow in ArcMap

but I have spent days, if not weeks digging through all the responses with basically little evidential results or improvements. For this reason I am asking again, and bringing to the fore.
When using ArcGIS Pro I experience an awful lot of slow speed issues, including:

time to take to click between ribbon tools,
opening tools,
running tools,
panning through maps with minimal or no data on,
folder operations,
cancelling geoprocessing tools.

Please note:

All data is local, except for the ESRI basemap. (it's faster without the ESRI basemaps but the software and ribbon use is still incredibly slow -  I don't even try to look at the Geography of Map layers in the Catalog view).
All data layers are transformed on the fly to the Projection system of the basemap (WGS)- this is only a short term solution as it shouldn't be a requirement; people need to to work in a lot of other projection and CRSs.
I have 60mbps bandwidth
My machine has 16GB of RAM and is and i7-9850 machine; its specified FOR GIS.

Through testing it has to be the software because when I open the same data in QGIS, it's almost instantaneous.
I don't really know what to do.


texts: conversion of *.tif file format to *.hdr file format

texts: I have clipped DEM of a region from SRTM data. It is in .tif format. Now I want the same file in header format (.hdr). How can I convert *.tif file to *.hdr file?


texts: Rasterize vector layers with different extents into identical raster grids with QGIS

texts:
I am working on a quite simple spatial potential analysis using QGIS. I have some ArcGIS knowledge, but I am new to QGIS.
My task is to identify the functional inner city of a region based on different criteria like density or distance to public transport. Therefore I created polygon layers representing different criteria. All the polygons have values between 0 and 1. The next step is to rasterize the vector layers and use the raster calculator to make an additions for each raster cell.

The problem I am facing is that the layer extents are different. Therefore all the raster outputs also differ in extent and cell size. As a result I cannot use the raster calculator to make the cell calculations.

So what I need to do, is to find a way to rasterize the polygone layers into raster sets with the same extend and the same cell size. ArcGIs offers the options to set a specific “Analysis extent” and a specific “Analysis cell size” for further raster calculations, which is basically what I am looking for in QGIS. I hope you guys can help me!


texts: R: calculating overlap between two standard deviational ellipses

texts:
Aim

Calculating the area of overlap of two standard deviational ellipses in R and work out their intersection and union.

Tool to be used

calc_sde() from the 'aspace' package: calc_sde()

Question

Given some of the values returned by the function (see below), is it possible to calculate the area of the ellipses and the extent of their overlap?
In particular, I would like to be able to work out the intersection and union between the two polygons. Not an expert in trigonometry and maths, so I am wondering which of the following values returned by the aboventioned function can be put to work and how:

CENTRE.x
X-coordinate of the centre

CENTRE.y
Y-coordinate of the centre

Sigma.x
Half-length of axis along x-axis

Sigma.y
Half-length of axis along y-axis

Theta
Rotation angle in degrees

Eccentricity
A measure of eccentricity (i.e., the flatness of the ellipse)

Area.sde
Area of the SDE

TanTheta
Trigonometric result

SinTheta
Trigonometric result

CosTheta
Trigonometric result

SinThetaCosTheta
Trigonometric result

Sin2Theta
Trigonometric result

Cos2Theta
Trigonometric result

ThetaCorr
Corrected theta angle for rotation of major axis from north


texts: shapefile not showing up in results after running a customized model

texts: I made a custom model containing two buffers and one erase to show the area within 2 miles, but 0.5 miles away from major roads. The model runs perfectly, but there isn't any result produced. The shapefile for both buffers are available in the catalog but the shapefile that should be produced by the final erase function is not there and could not be found in the folder either. Here is my model.


texts:
I was hoping to use ESRI's ArcGIS Editor for OpenStreetMap to view OpenStreetMap data in ArcMap, but my organization hasn't upgraded to version 10 yet.

I also know about ArcBruTile, but I prefer not to use that software because of the potential TOS violations with Google, Bing, etc.

Does anyone know of another way to bring OpenStreetMap tiles into ArcGIS 9.3.1?


texts:
I've discovered many annoyances with Leaflet.js related to rendering markers. They frequently cover each other, even when zoomed in fully, forcing me to balance the cursor with 1-pixel precision to get to the pop-up data (or access for clicking) the marker that is "behind" another marker. It seems to place itself 1 pixel apart to be "helpful", but this is not nearly enough.

Is there some way to get any markers that overlap to start "bending" to the side, so that they still point to the right position, but still are "out of the way" so you can actually easily click them (not to mention SPOT them!)?

I assume that this feature already exists, or something like it, because it's one of the first things I've noticed having dealt with maps for not very long.


texts:
I don't know why the users that visit my website cannot see the map i have published with geoserver.
The steps i took to do this have been:

importing the data to geoserver
adding the layers to geoexplorer
i have changed the style a little
hit publish map on the geoexplorer
selected the tools i want to have - next
copied the link that i need to copy in the webpage
changed from localhost to the name of my pc
pasted the code in my webpage.
i can see the page with the map frame, but the visitors cannot.
the error they get is:

  the server refused the connection


texts: how to convert latitude, longitude projection from epsg:3857 to epsg:4326

texts:
I want to use geocode in our application. so that, I got Google geocoding API and could achieve with Google Maps.
Problem is that I wanted to achieve the same with my own map which is having projection EPSG:4326 INSTEAD OF getting geocoded from Google Map which is having projection EPSG:3857.
or
 I am getting a coordinates from Google geocoding projection EPSG:3857 and that
coordinates I want to mark on my OpenLayers Map which is having projection EPSG:4326

please help me.


texts:
in the  'System Variables' i set

Path =C:\OSGeo4W64\apps\qgis\bin

This is where the 'qgis_core.dll' is located.

PYTHONPATH=C:\OSGeo4W64\apps\qgis\python

But when i import qgis.core i get:

   run o-help for a list of available commands
C:\>python
Python 2.7.4 (default, Apr  6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> import qgis.core
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\OSGeo4W64\apps\qgis\python\qgis\__init__.py", line 27, in <module>
    from qgis.core import QgsFeature, QgsGeometry
  File "C:\OSGeo4W64\apps\qgis\python\qgis\core\__init__.py", line 30, in <module>
    from qgis._core import *
ImportError: DLL load failed: %1 is not a valid Win32 application.
>>>

When I change the paths to osGeo4W directories (Rather then the osGeo4W64) I get:

run o-help for a list of available commands
C:\>python
Python 2.7.4 (default, Apr  6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> import qgis.core
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\OSGeo4W\apps\qgis\python\qgis\__init__.py", line 27, in <module>
    from qgis.core import QgsFeature, QgsGeometry
  File "C:\OSGeo4W\apps\qgis\python\qgis\core\__init__.py", line 30, in <module>
    from qgis._core import *
ImportError: DLL load failed: The specified module could not be found.
>>>

Can sombody understand what I'm doing wrong?


texts: How to use Select by attributes more efficiently

texts:
I made a shapefile into 1,040,000 grids by the tool, create fishnet. This shapefile has an attribute table including 7 fields and 1,040,000 rows. What I want to do is to make a model which can let me get some "values" that I would like to extract.( I am sorry that I don't know how to describe these "values" in a sentence. I will explain it in the following description.)

Its fields are like this:

Aspect (Long): (with values) 1, 2, 3, 4, 5

Slope (Long): (with values) 1, 2, 3, 4, 5, 6, 7, 8, 9, 10

Curvature (Long): (with values) 1, 2, 3, 4, 5, 6, 7

Rainfall 1hr (Double): (with values) 89.9, 100, 120, 140, 160, 160.1

Geology (Text): (with values) RockA, RockB, RockC, RockD, RockE

Distance (Long): (with values) 0, 1000, 2000, 3000, 4000, 5000

Rainfall 3hr (Double): (with values) 239.9, 320, 400, 480, 560, 640, 720, 800, 800.1

I am trying to fill in the tables below(this picture does not include all the tables that I want to fill in. It is just part of all the tables):

For the yellow cell in the first table. I should fill in the row number,or you can say the grid number, with attribute "Aspect = 1" and "Slope = 1". For other cells in those tables, they have the same meaning.

I will do this search for all the 7 fields two by two. So, I will need to complete [(6+1)*6]/2 = 21 tables.

Up to now, I have built a simple model that uses tool Select Layer by Attribute, Get Count, and Collect Values. By using this simple model, I can get the values of a column in one search. But I think it is still too time-consuming.

The below pictures indicates the values extract from one of my model. After running my model I will type these numbers into my table by hands.

I am also curious if one can extract these values into .txt or excel files. Then, it will be much easier and convenient for me. I know there is a similar question here(Exporting data from Collect Values output in ArcGIS ModelBuilder?). I have tried it successfully. However, I can only get the values for "one column" written for each search, similarly.

I hope someone could teach me how to do this search more efficiently, such as getting the values for two fields in one search.

I hope my explanation is clear enough. Any suggestions are welcome.


texts:
I have a feature set of locations, made up of several points.

I have the Network Analyst extension for ArcMap, and a feature set of all the streets in the area I'm looking at.

How can I create drive time polygons (polygons of the areas I can reach, for each point feature, within a certain number of minutes) using these tools?

So far I've only found tutorials for doing this with the Business Analyst extension, which I do not have and is not available on MyESRI.


texts: Do we have a way to change the colors of the polygons directly in JOSM?

texts: Do we have a way to change directly the colors of the polygons in JOSM?


texts: Accuracy Assessment for single date classification - possible to apply to classification without reference images?

texts:
I applied supervised maximum likelihood classification to a yearly stack of Landsat images (resolution 30m). The land cover classes are related to coniferous forest, i.e., forest, clear-cut, fire, bark beetle. I had only few years of available reference aerial photos (resolution 0.5m).
Aerial photos do not fully cover the extent of Landsat images (AOI). To assess the accuracy of the classification for single year scene, I extracted a overlapping part between the extent of Landsat and aerial photos:

I calculated the measures of accuracy for the extracted part of fully classified Landsat by set of stratified sampling points per mapped class.

I followed the methods published here:

Olofsson, P., Foody, G.M., Herold, M., Stehman, S. V., Woodcock, C.E., Wulder, M.A., 2014. Good practices for estimating area and assessing accuracy of land change. Remote Sens. Environ. 148, 42–57. doi:10.1016/j.rse.2014.02.015 and
Olofsson, P., Foody, G.M., Stehman, S. V., Woodcock, C.E., 2013. Making better use of accuracy data in land change studies: Estimating accuracy and area and quantifying uncertainty using stratified estimation. Remote Sens. Environ. 129, 122–131. doi:10.1016/j.rse.2012.10.031

to assess the accuracy of the classification.

Great R code summarizing this can be found here:
https://github.com/openforis/accuracy-assessment/blob/master/error_matrix_analysis.R

The error matrix, by this method, is not based on the simple sample counts, but on the estimated proportions of the area. Including the proportion of the mapped classes in accuracy assessment, it is possible to quantify the uncertainty attributable to variability in the sampling. Simply, to adjust the mapped area to map classification error. The error-adjusted area estimates should be accompagnied by confidence interval.

My problem:

I classified 7 Landsat images. I calculated accuracy for extracted parts of them for two years, because of no availability of another reference data sets.

My questions:

How to calculate the adjusted mapped area for the whole classified Landsat scene? (i.e. the overlapped area, where I assessed accuracy, was 100 ha. The whole AOI was 500 ha. My adjusted area estimated for the overlapped area is 20 ha (100ha±20 ha)).
How to extent the margin of error number for the whole classified area? (500±100 ha?)