texts: Automated way to check if layers are exactly the same using QGIS?

texts:
I'm doing a voluntary project where I'm helping a client with an objection to a new road construction. The client, without much GIS knowhow, has provided me with a load of .dbf, .prj, .sbn, .sbx, .shp, .gfs and .shx-files he dragged from the municipality website about the project.

From a first visual inspection, the files seem to contain the same data, but I want to be sure before tossing them out; is there any way to have QGIScheck if the layers are the same?


texts: Reducing Resolution using Google Earth Engine?

texts:
I am just starting to use GEE to do some raster preparation for ecological niche modeling. Mostly, I am using MODIS images at 1,000 m resolution. This is the scale I need all my rasters in. MODIS uses sinusoidal projection "SR-ORG:6974". Now, in processing SRTM Digital Elevation Data Version 4, I performed raster export after selecting my export area as follows:

// Display DEM model
var DEM = ee.Image('CGIAR/SRTM90_V4')
.select('elevation')
.clip(Argentina)
print('DEM', DEM.projection().nominalScale())
print('DEM', DEM)

Map.addLayer(DEM, {}, 'DEM')

// Create export region
var ExportArea = ee.Geometry.Rectangle([-73.7,-55.0,-53.4,-21.4]);
Map.addLayer(ExportArea, {color: 'FF0000'}, 'poly');

// Export raster 1
Export.image.toDrive({
image: DEM,
description: 'DEM',
folder: 'My drive folder',
crs: 'SR-ORG:6974',
maxPixels: 1e13,
scale: 1000,
region: ExportArea,});

I need a final raster at a resolution of 1,000 m. Although the above procedure is straightforward, I was wondering if it is the best approach to reduce resolution from 90 (DEM) to 1000 m (MODIS).

Then, following GEE developers' advice on how to aggregate pixels to larger pixels (MODIS at 1,000 m) in a different projection https://developers.google.com/earth-engine/resample I did the following:

// Load a MODIS EVI image.
var modis = ee.Image(ee.ImageCollection('MODIS/006/MOD13A1').first())
    .select('EVI');

Map.addLayer(modis, {}, 'MODIS')

// Get information about the MODIS projection.
var modisProjection = modis.projection();
print('MODIS projection:', modisProjection);

// Load and display DEM data at 90 m resolution
var DEM = ee.Image('CGIAR/SRTM90_V4')
.select('elevation')

print('DEM scale', DEM.projection().nominalScale())

// Get information about the DEM projection
var DEMProjection = DEM.projection();
print('DEM projection:', DEMProjection);

// Get the DEM data at MODIS scale and projection
var DEM_reproj = DEM

// Force the next reprojection to aggregate instead of resampling
.reduceResolution({
reducer: ee.Reducer.mean(),
maxPixels: 1024
})

// Request the data at the scale and projection of the MODIS image
.reproject({
crs: modisProjection
});

// Create region
var ExportArea = ee.Geometry.Rectangle([-73.7,-55.0,-53.4,-21.4]);
Map.addLayer(ExportArea, {color: 'FF0000'}, 'poly');

// Export raster 2
  Export.image.toDrive({
  image: DEM_reproj,
  description: 'DEM',
  folder: 'WildBoar',
  maxPixels: 1e13,
  scale:1000,
  region: ExportArea,});

After exporting both of the above raster to My Drive, I explored them in QGIS and obtained the following information for each:

Exported raster 1 (QGIS)

Extent
-73.7067690620067424,-55.4260530301744723 :
-53.3958604880643577,-21.3978700677270020

Unit: degrees

Origin: -73.7068, -21.3979

Pixel Size: 0.00898315, -0.00898315

Exported raster 2 (QGIS)

Extent
-7631000.0000000000000000,-6163000.0000000000000000 :
-3405000.0000000000000000,-2379000.0000000000000000

Unit: meters

Origin: -7.631e+06, -2.379e+06

Pixel Size: 1000, -1000

So at this point my main question is: Which method is the correct one? As seen in the screenshots above, the first script outputs a raster with pixel size 0.00898315 by 0.00898315 degrees, which equals 997.12965 m E/W at 45° S (center of Argentina). Exproted raster 2, on the other hand, consists of a raster with a 1,000 by 1,000 m pixel size. The way pixels are aligned in Output 2, however, looks odd to me, as evidence by their inclination. I am not sure what may account for this, or if this is an expected behavior at all. By the way, pixels appear tilted in GEE as well.


texts: What does "faceting" mean in a geo/map context?

texts: I saw some folks talking about "faceting" in relation to a map data display, what does this mean in general terms? Examples?


texts: Find Properties of Nearest Line Segments to Point in Python

texts:
I have been trying to find a workable Python approach to finding nearest line segments to one (or more) points, and came accross @gene's solution to @RustyMagoo's problem here:

from shapely.geometry import Point, LineString, mapping, shape
from shapely.ops import cascaded_union
from shapely import affinity
import fiona

sweep_res = 10  # sweep resolution (degrees)
focal_pt = Point(0, 0)  # radial sweep centre point
sweep_radius = 100.0  # sweep radius

# create the radial sweep lines
line = LineString([(focal_pt.x,focal_pt.y), \
                   (focal_pt.x, focal_pt.y + sweep_radius)])

sweep_lines = [affinity.rotate(line, i, (focal_pt.x, focal_pt.y)) \
               for i in range(0, 360, sweep_res)]

radial_sweep = cascaded_union(sweep_lines)

# load the input lines and combine them into one geometry
input_lines = fiona.open("input_lines.shp")
input_shapes = [shape(f['geometry']) for f in input_lines]
all_input_lines = cascaded_union(input_shapes)

perimeter = []
# traverse each radial sweep line and check for intersection with input lines
for radial_line in radial_sweep:
    inter = radial_line.intersection(all_input_lines)

    if inter.type == "MultiPoint":
       # radial line intersects at multiple points
       inter_dict = {}
       for inter_pt in inter:
           inter_dict[focal_pt.distance(inter_pt)] = inter_pt
       # save the nearest intersected point to the sweep centre point
       perimeter.append(inter_dict[min(inter_dict.keys())])

    if inter.type == "Point":
       # radial line intersects at one point only
       perimeter.append(inter)

    if inter.type == "GeometryCollection":
       # radial line doesn't intersect, so skip
       pass

# combine the nearest perimeter points into one geometry
solution = cascaded_union(perimeter)

# save the perimeter geometry
schema = {'geometry': 'MultiPoint', 'properties': {'test': 'int'}}
with fiona.open('perimeter.shp', 'w', 'ESRI Shapefile', schema) as e:
     e.write({'geometry':mapping(solution), 'properties':{'test':1}})

What I need in addition is to extract the nearest lines' attributes/ properties, like id, name, shape, etc. For example, if my shapefile of lines are rivers, and attributes of the rivers are name, length, data about the amount of water passing each day, etc. How would I go about extracting this information using the method proposed by @gene/@RustyMagoo?


texts: Google Earth Engine compute point-polygon distances

texts:
I am trying to compute distances from points stored in a fusion table to the nearest political border (also from a Table) in Google Earth Engine.

I need to compute distances between points.  I would need the country name at the nearest border (it is stored under the variable 'name' in the border-table).

Please find the code below,

var fg_points = ee.FeatureCollection('ft:1DVI9dvjC00QhP-
hZgO3wQeH45RO2FDtdVcvF0fZq');
var borders = ee.FeatureCollection('USDOS/LSIB/2013');
Map.setCenter(74.582748, 42.882004, 8);

// Define a spatial filter with max distance to points is 500.
var spatialFilter = ee.Filter.withinDistance({
distance: 500000,
leftField: '.geo',
rightField: '.geo',
maxError: 10
})

//join the points to borders
var joined = ee.Join.saveAll({
matchesKey: 'neighbors',
measureKey: 'distance',
ordering: 'distance'
}).apply({
primary: fg_points,
secondary: borders,
condition: spatialFilter
});

// Get rid of points w/o neighbors.
var hasNearest = joined.map(function(f) {
var neighsSize = ee.List(f.get('neighbors')).size();
return f.set('neighsSize', neighsSize);
}).filter(ee.Filter.gt('neighsSize', 1));
Map.addLayer(hasNearest, {color: 'red'}, 'hasNearest');

// Get distance to nearest point.
var withNearestDist = hasNearest.map(function(f) {
var nearestDist = ee.Feature(ee.List(f.get('neighbors')).get(1))
.get('distance');
return f.set('nearestDist', nearestDist);
});

Export.table.toDrive(hasNearest,
"Distance",
"Distance",
"Distance");


texts: Using exponential integral function for raster calculation in ArcPy

texts:
How do I use "exponential integral" function for raster calculation in ArcPy?

The equation that I want to apply in raster calculator is of the form EXPINT (r^2*S/(4*T*t)). Here, r, S, and t are constant parameters and I have raster file (say raster1 ) for T parameter.

I am a beginner.

import arcpy
import numpy as np
from scipy import special as sp
from arcpy import env
from arcpy.sa import *

raster1 = ("D:/raster1")
raster2 = ("D:/raster2")
raster3 = (raster1 * raster2 * 30 * 4 * np.pi)/-sp.expi(- (0.20 ** 2 * 0.0001/(4 * raster1 * 30)))


texts: In which projection I should project a polygon of a country to get the right area in the attribute table without a distortion? For example greenland


texts: Sharing custom made color ramps in ArcGIS Desktop?

texts:
I've successfully created a custom color ramp in ArcMap (v. 10.3). From the Style Manager I can see the path to where my custom color ramp is located, but I'm unable to locate that path on my machine; it's like it doesn't exist. I would like to be able to share the color ramp I created with other coworkers, but can't figure out how to locate the path to the color ramp or how to share it.

How can I resolve this issue?


texts: Merging SHP lines with attribute data combined?

texts:
I have one shapefile which contains hiking routes and another shapefile which has ski routes. Sometimes they follow the exact same path.

In QGIS if the lines mentioned above are in the same SHP I can merge two lines in edit mode by selecting them and Merge selected features and then chose Concatenation with an attribute field so that the value of the merged line will be "hiking,skiing" etc.

This method takes time since I have thousands of lines like this.

Is there a way to do this for whole shape files? Maybe five or six of them.

I saw How can i merge lines with same name and intersect each other? which is similar but I'm not familiar with spatialite and would like to stick to GDAL and/or QGIS.


texts: Calling Gdal contour from python/ipython

texts:
Hello I am trying to patch process a folder full of ASCs to contours using Gdal from python. here is my code;
import os
import glob

PATH1 = raw_input("folder containing ASCs path please:")

inputlistB = []

for filename in glob.glob(os.path.join(PATH1, '*h_Max.ASC')):
    print filename
    inputlistB.append(filename)

print inputlistB

for filepath in inputlistB:
    src = filepath
    dst = filepath+'b.shp'
    os.system('gdalcontour -fl -888 -inodata -f "ESRI Shapefile "' + src + " " + dst)

The ipython console returns this;

folder containing ASCs path please:C:\D\OUT
C:\D\OUT\Wray_0010F_012_h_Max.asc C:\D\OUT\Wray_0020F_012_h_Max.asc
C:\D\OUT\Wray_0030F_012_h_Max.asc
['C:\D\OUT\Wray_0010F_012_h_Max.asc',
'C:\D\OUT\Wray_0020F_012_h_Max.asc',
'C:\D\OUT\Wray_0030F_012_h_Max.asc']

showing that the inputlist functions as intended. However no new gis layers are generated in my target folder. If I take the statement ran by os.system and copy it with file paths into command line, it works. Though it does crash after producing the shapefile. If I just run the final line of the for loop os.system(... the Ipython console simply returns :1
Any idea why gdal would crash after producing a shapefile / successfully running?
Any ideas as to why my code does not produce any results despite no python errors?
Any alternative / better solutions to this problem?
also I use the file path input and loop in a lot of programs, feel free to steal it if this is news to you.


texts: Repository for GPS tracklog (better gpx format)

texts:
I've developed some algorithms that work with GPS trajectories. So I need a repository in which I can download GPS trajecotries to run my test.

I most need marine and air tajectories, better if fileformat was a .gpx.


texts: QGIS: How to modify the selection style when selecting features in a standalone script

texts:
It's quite easy to modify the default style that QGIS applies when selecting attributes in the GUI by following the steps in this question: QGIS: How to modify the appearence of selected attributes
However, when working with a stand-alone script that opens the same project, selects a feature and exports the contents of the bounding box to a png, the selection style is still the same default yellow.
For rendering and saving, I'm following the steps described in the QGIS Cookbook. The only addition is feature selection and setting extents to feature bounding box.
I probably have to set the style explicitly in a stand-alone script. How to accomplish that?
EDIT: MINIMAL REPRODUCE
Step 1: modify selected feature style in project
Step 2: load project and try to render layer with the following code
project = QgsProject.instance()
project.read('/where/my/project/at')

vlayer = project.mapLayersByName('MY_LAYER')[0]
vlayer.selectByIds([my_id])
settings = QgsMapSettings()
settings.setLayers([vlayer])
settings.setBackgroundColor(QColor(255, 255, 255))
settings.setOutputSize(QSize(800, 600))
settings.setExtent(vlayer.extent())

render = QgsMapRendererParallelJob(settings)
    image_location = os.path.join(QgsProject.instance().homePath(), "render.png")

def finished():
    img = render.renderedImage()
    # save the image; e.g. img.save("/Users/myuser/render.png","png")
    img.save(image_location, "png")

render.finished.connect(finished)

# Start the rendering
render.start()

# The following loop is not normally required, we
# are using it here because this is a standalone example.
from qgis.PyQt.QtCore import QEventLoop
loop = QEventLoop()
render.finished.connect(loop.quit)
loop.exec_()

Step 3: feature selected but selection color still standard QGIS yellow.


texts: Satellite image in ArcGIS darker than usual

texts:
This is one image from 2013 in some área.

This is another image from 2014 in the same área, but darker. why? How can I display it like in 2013.

Both are Landsat 8 images.


texts: Reader character encoding problems in Oracle layers with where parameters in FME Desktop

texts:
I am using FME Desktop 2020 and trying to read an Oracle layer with -where- condition but i have problem because it failed with reading of character Ë or Ç.

I get this error:

ORA-01756: quoted string not properly terminated'

If I replace Ë with E, it is ok. Find below part of the error log.


texts: Setting list of fields to field mapping of spatial join using ArcPy?

texts:
I try to code spatial join tool and delete some fields in the result. to delete the fields i list the fields of source layer and the result layer then use the list of differences in the two list and  usedelete field management tool.it's working but not good solution.If i can set the Spatial Join field mapping to remain the source layer fields and only the "Distance_To_Next" field of the join layer, i get rid of some codes and iterates but don't know how to set the field mapping in the code.

Scripting spatial join in ArcPy to set field attribute in field mapping object? is the only result that i could find in the stack but not my answer.

import arcpy
arcpy.SpatialJoin_analysis("Point1", "point2","JoinResult","JOIN_ONE_TO_ONE","#","#","INTERSECT")
listfields1 = [fraft.name for fraft in arcpy.ListFields("Point1")]
listfields2 = [fraftv.name for fraftv in arcpy.ListFields("JoinResult")]
diff_fields = [i for i in listfields2 if not i in listfields1]
for field in diff_fields:
    if field != "Distance_To_Next":
        arcpy.DeleteField_management("JoinResult",field)


texts:
So I'm trying to use the quiet option for gdal calc, this is a cut of my code:

    call([sys.executable, 'C:\Program Files (x86)\GDAL\gdal_calc.py', '-A',
    Outer_Output_Raster, "--outfile=" + New_Outer_Raster, '--calc=-9998*
    (A>-9999)','--NoDataValue=-9999','--overwrite','--quiet'])

And I get the following error:

    gdal_calc.py: error: no such option: --quiet

But the documentation clearly lists this as an option, see documentation here.

Does anyone have any suggestions or am I missing something completely?

For the record the rest of the code works exactly as intended, it is just the quiet option that isn't working.


texts: PyQGIS: drop features belonging to same group conditional on other attribute expression

texts:
I have a polygon shapefile on which I need to eliminate features belonging to the same group (identified by a variable/attribute XG) on the base of the values taken by another variable YA. In particular, given some features with same XG, I would like to keep the one with the highest YA value.
There are many groups for which XG identifies only one feature.

Below is a sample structure of my attribute table in case I did not explain my self well.
What I have:

XG YA
5 100
5 130
5 210
2 0.4
2 0.1
1 0.6

What I would like to get using Python console in QGIS:

XG YA
5 210
2 0.4
1 0.6


texts:
I have a city wide water distribution network that will undergo a flushing process over the next couple of years. I need to keep track of what segments of pipe have been flushed. The work is divided into zones, and I know what pipe unique ID is in each zone.

Say that I have one segment of pipe that belongs to multiple flushing zones. How could I design a field that allows one attribute to belong to multiple zones, so that when I indicate if a zone is complete it gets updated with the most recent completion date?

Are there any ways to make a field like this? Then I am thinking I could write  a query that says IF: (unitID) belongs to zone (1,2,3, etc.) THEN update new field with the given date.

Any thoughts on this?


texts: Calculating Geometry (area) in a spatial join attribute table

texts:
In ArcMap, will Calculate Geometry (area) be accurate for records that have other polygons attached via Spatial Join?  It doesn't seem to be adding the extra polygons that are spatially joined. ArcGIS 10.

However, my original Target shapefile, USFS Land Status layer, might have not excluded the areas for Fresh Water lakes in the first place even though the original attribute table breaks them out separately.

So either Spatial Join isn't adding the Fresh Water lakes in or the original Land Status Layer is counting the areas twice 1) within the area calculated for the land records and 2) in the Fresh water records.  I had broken the original records apart via the Select tool and then did Spatial Join.


texts:
I have a simple question, but I can't find the answer
I am looking for the principal angle of a minimum_rotated_rectangle polygon in relation to a latitude or a longitude.
df4.minimum_rotated_rectangle

Does anyone have this in stock?


texts: Attempting to run first script after setting up environment for QGIS and python?

texts:
I uninstalled arcmap and all existing python installations. Then I followed the instructions from https://www.youtube.com/watch?v=Chq9fTRKsHk&t=2s to install qgis and set the environment variables.

I'm attempting to run the script from the second video (https://www.youtube.com/watch?v=Y8JjvBKOCDk&list=PLzHdTn7Pdxs4AvkJeJwlyav2-at9QwKWR) which should launch a blank canvas.

from qgis.core import *
from qgis.gui import *
from PyQt4.QtCore import *
from PyQt4 import QtGui

app = QtGui.QApplication

QgsApplication.setPrefixPath("C:/OSGeo4W/apps/qgis", True)
QgsApplication.initQgis()

window = QtGui.QMainWindow()
window_frame = QtGui.QFrame(window) 11
window.setCentralWidget(window_frame)
frame_layout = QtGui.QGridLayout(window_frame)

canvas = QgsMapCanvas()
frame_layout.addWidget(canvas)

window.show()
app.exec_()

It's supposed to bring up a canvas but does not.

The error I'm getting is

Process finished with exit code -1073741819 (0xC0000005)

I don't know if the problem is in the environment variables or in the script. Any ideas?


texts: Selecting two or more values (like two city names) in one column

texts:
In the attribute table I have a column with many cities and want to select rows that only contain these city names, like "Paris" and "London".
I tried: "cityname" ILIKE 'Paris' OR 'London' but it doesn't seem to work.
I have seen a reply using "arrays" to a similar question but I am not sure how it works and I would prefer to use ilike if possible.
Selecting more features with the "ilike" and "in" operator using QGIS


texts: netCDF does not display properly in QGIS 1.8

texts: I have a netCDF file of fire emissions on a regular grid of 12 km by 12 km grid cells with 200 east-west grid cells by 161 north-south grid cells centered on the continental United States.  The map projection is Lambert Conformal.  I am using QGIS v1.8.  When I open the raster it displays as a very, very tiny raster located in mid-northern reaches of the State of Kansas.  Can someone tell me what blunder I am performing?


texts: Is PDF version exported from ArcGIS Pro different to PDF from ArcMap?

texts:
Is there a difference between the PDF version exported from ArcGIS Pro 1.4.1 layouts and that produced by ArcGIS Desktop 10.4.1?

I ask because I am seeing PDFs exported from ArcGIS Pro (600dpi) that look perfect in Adobe Acrobat Reader DC but when sent to either of two printers (no others tested) there are some grid lines (from Create Fishnet) that show up about 10 times as thick as they should.  The same thing happens in the scale bar and there is bad "blockiness" in two picture elements (logos) on the layout too.

Although the two PDFs are created quite differently (ArcGIS Pro vs ArcMap, using ArcPy in each case to do the same map series) these symptoms are not seen in the ArcMap PDFs but are seen in all such maps from Pro.

The hardcopy is from an A0 HP plotter and an A3 printer so I think the problem is the PDF rather than either plotter/printer.

A mobile phone photo of the hardcopy, although quite blurry, illustrates the symptoms.  Note the thick black line on the scale bar at bottom right (the rest is the expected thin lines), the two thick blue lines on the index map at lower left (the rest are the expected thin lines) and two thick grey lines on the grid in the main map (the rest are the expected thin lines).

Curiously, the same layout exported as a TIFF (600dpi) and imported into Adobe Acrobat Pro before being saved as a PDF prints fine from that PDF.  However, this is a workaround that is only viable for a subset of requirements because it requires manual creation of the final PDF (in an otherwise automated system), results if TIFF file sizes that are 10x as big and take 10x as long to create, and require licensing of Adobe Acrobat Pro.


texts: OSM2po is converting data from OSM to postgresql. But how is it handled with places or squares? Not every places, square has pedestrian way routes on it. But for routing I need nodes and pedestrians can walk through places without restriction, if there are no obstacles between. But we can consider only "empty" places. Therefore the shortest distance is to walk through the place. How is osm2po handling area=yes or highway=footway for areas, respectively.


texts: Is it possible in QGIS 2.18 to change the order of the pages in the composer (eg. move the last page to page 4)?


texts: Extracting specific roads from raster image using ArcGIS Desktop?

texts:
I want to create a line shapefiles of white roads from the raster file (screenshot1) using ArcGIS:. My idea is to use ArcScan toolbar, but first of all I need to classify roads pixels values from raster. For that I have used ISO Unsupervised classification tool with a number of 30 classes. The results are not good for me because I got the unclear pixel values which are unusable for ArcScan vectorisation (screenshot2):

Of course, if I have a clear pixels defining the roads I will use Reclassify tool to filter get just two pixel values (one which I need and one which don`t).
Now I need to know is there any way to solve my problem?


texts: is possible data management with geoserver rest api

texts: We can create data stores, workspaces with goserver rest api. Is possible to load data to a datastore via geoserver rest api?


texts: Opening Attribute Table by using ArcGIS Pro SDK, WPF and C#

texts:
I am trying to display an attribute table for a selected layer. I have access to attribute table but I can not pass the ArcGIS table to the DataGrid I created in xaml.

<DataGrid x:Name="DataGrid_A_Table" HorizontalAlignment="Left" Height="100"
Margin="10,48,0,0" Grid.Row="1" VerticalAlignment="Top" Width="280"/>

Do you have any idea?

public void create_Layer()
        {
            string shp_path = this.ComboBox_shpfiles.Text;
            int indexNumber = 0;

            if (shp_path != null && shp_path != string.Empty)
            {
                System.Uri shp_filepath = new System.Uri(@shp_path);
                QueuedTask.Run(() =>
                {
                    var shp_layer = LayerFactory.Instance.CreateLayer(shp_filepath, MapView.Active.Map, indexNumber, shp_path.Substring(16));
                    //Get the active map view.
                    var mapView = MapView.Active;
                    if (mapView == null)
                        return Task.FromResult(false);
                    else
                    {

                        var selectedLayer = mapView.GetSelectedLayers()[0];
                        //var tableContent = TableControlContentFactory.Create(selectedLayer);

                        if (selectedLayer is ArcGIS.Desktop.Mapping.FeatureLayer)
                        {
                            ArcGIS.Core.Data.Table a_table = (selectedLayer as FeatureLayer).GetTable();

                            DataGrid_A_Table.DataContext = a_table; // This line is not working

                            int gugu = a_table.GetCount(); //I can the count of rows here
                            MessageBox.Show(gugu.ToString());
                        }

                        //Zoom to the selected layers in the TOC
                        return mapView.ZoomToAsync(selectedLayer);
                    }
                });
            }
        }

I tried to cast ArcGIS.Core.Data.Table to System.Windows.Documents.Table but did not work

I am getting wrong thread exception.


texts: How to run "create mask" from the console?

texts:
There's a plugin called "create mask" for qGIS and trying out some simple scripting from the Python console.

I would like to call it after I select the vector features I want, then do more related tasks afterwards. So, how does one run "create mask" from the console?


texts: How to interpret Orfeo-Toolbox Hyperspectral Unmixing results and calculate percentage of endmember abundance?

texts:
I understand that using the Hyperspectral Unmixing Tool in the Orfeo-Toolbox gives me the abundance matrix of the endmembers in a given pixel, what I want to do is to calculate the percentage each endmember contributes to that pixel.
An example pixel from my work would be: [137.75, 123.70, 58.39] which results in the following abundance matrix with 4 endmembers: [0.76, 0.03, -0.35, 0.92].
As far as I understand, these are the coefficient of the linear equation, that when multiplied with the end members gives me the pixel value. How can I now calculate the percentage of the abundance each endmember has in the pixel?


texts: Geographically Weighted Regression and categorical values leaving blank values after analysis

texts:
I am using Geographically Weighted Regression (GWR) in ArcMap 10.2 to look at the effect of environmental and human factors on lion occurrence in Namibia.

I created a vector grid over the study area, calculated the frequency of occurrence of lion GPS points, calculated the distance to each nearest feature using the near tool from the analysis toolbox and "coded" each grid cell with the categorical values.

When I run a GWR to look at the effect of say distance to river it works fine and I have found the suitable bandwidth for the analysis.

When I run a GWR to understand the effect of land use (i.e. national park or not national park) it calculates the values along the edges of the parks boundaries but the grid cells in the park come out blank after the GWR run and in the attribute tables it says null. I coded each vector grid cell with "1" if it is inside the national park and "0" when it was outside. When I input the data into GWR, the dependent variable was the lion frequency of occurrence and the explanatory variable is landuse (which is the the 1 and 0's for park and non park). I also used the near tool to calculate the distance to the park boundary but this gave similar results when using GWR and gives everything inside the park a value of 0.

The lions do spend most of their time in the park and that is why I wanted to show that the fact that the park there is a contributing factor that lions still occur in that area.

Why would the values come out blank? If it helps to explain the situation better I can also upload an image of the results.


texts: Creating a contour raster from contour lines shapefile

texts:
I'm trying to create a 3D view in Qgis 3.0 for the first time. It's asking for a raster file from which to calculate elevations. All I have for elevation data is a shapefile of contours (as lines).

I've tried simply rasterizing the contour lines shapefile, but it comes up with Error 1 "Wrong value for -outsize parameter" and I can't seem to get any information on this error.

Is there a way to get a raster file suitable for use with the 3D view from a shapefile with contour lines?

Edit - I've followed the link trail from the previous question, they all seem to link back to the tutorial which @Marco linked to.

The trouble is, all those start with converting the vectors to raster, and that, as the question suggests, is where I'm having the problem. Perhaps, given that converting raster contour lines to raster DEM is so well covered I should change the question title to just "How do I convert vector contours to raster contours?"


texts: GDAL Slope function not working with SRTM DEM data?

texts:
I'm a QGIS rookie trying to generate slope data for research projects in South Africa.
I've managed to download DEM data from Bing Aerial maps using the SRTM-downloader but cant get the GDAL slope function output in QGIS to work as seen in other successful efforts online. The map generated looks like the attached pic.The slope value ranges between 0 (black) and 89.99 (white).
The aspect function seems to work but can't crack the slope function, any suggestions?


texts: PostGIS algorithm to unite points of two geometries that are within specified radius

texts:
I want to find an PostGIS algorithm to unite points of two geometries that are "near" to each other.

That is: I want a function, that takes as input two geometries, and returns 1 geometry, that consist of points of first geometry, second geometry, and points that are in specified radius of both geometries.

Currently I am using such function for this (tolerance radius is 2 meters):

SELECT
ST_AsText(
ST_Buffer(
ST_Union(
ST_Buffer(
ST_GeomFromText(
'POLYGON ((0 0, 0 4.5, -2 4.5, -2 5.5, 0 5.5, 0 10, 1 10, 1 0, 0 0))', 3301
), 2
),
ST_Buffer(
ST_GeomFromText(
'POLYGON ((2 5, 5 8, 5 2, 2 5))', 3301
), 2
)
)
, -2)
);

It works reasonably well to unite points of two geometries, that are near to each other.

However, it has side effects. For example:

Input geometries:

Resulting geometry:

As can be seen from picture, some extra area is added to the "blue" geometry where the obtuse angles are. This is so, because using

ST_Buffer(ST_Buffer(geom, tolerance), -tolerance)

adds extra space to places where there are obtuse angles.

Is there an algorithm that would produce my asked result? Exact nature of how the two geometries are "glued" to each other is not so important. Buffering algorithm leaves smooth curve, but something like this would also be valid result (hand-drawn, I tried to add points that are in distance of 2 meters from both geometries):

If there are many algorithms, the fastest one will best suit my needs. But for comparison, slower algorithms will also be welcome (up-voted), if they excel by some other criteria.


texts: I need to include in the file containing geographical coordinates of my region a new external variable that I would like to use to grade the colour of the different parts of the whole map. I know how to include a new column but the only way I found to copy in my new variable is to write data by hands. Is that the only way?


texts: Add a common property to all sub-features of a feature collection in Earth Engine

texts:
I am new to earth engine and am struggling with a very simple task. I want to calculate the average NDVI for the counties of Maine for the month of January 2018. I manage to do so but, before exporting it to .csv, I wanted to add a property to each of the sub-features (counties) with the date ('2018-01-01') of the NDVI.

I try it with the .set() and setMulti functions but do not manage to achieve what I want. In particular, these functions add a property to the "overall" feature collection, and not to each sub-feature of the former. When I do that and export it as a .csv, the information is lost.

In short: I would like to add a property {date: '2018-01-01'} to each of the sub-features. Namely that the property in figure 1 is actually inside the properties in figure 2. I am attaching the code after the images.

Figure 1:

Figure 2:

Code:

// Load a collection of US counties with census data properties.
var maine = ee.FeatureCollection('ft:1S4EB6319wWW2sWQDPhDvmSBIVrD3iEmCLYB7nMM');

// let us use only counties from Maine:
maine = maine.filter(ee.Filter.eq('StateName', 'Maine'));

// loading recent NDVI data:
var dataset = ee.ImageCollection('MODIS/006/MOD13Q1')
                  .filter(ee.Filter.date('2018-01-01', '2018-01-31'));
var ndvi = dataset.select('NDVI');
ndvi = ndvi.mean();  // monthly mean -- need to average out as there are 2 images in this time period

// clipping into my AOI:
var ndvi_maine = ndvi.clip(maine);

// Average NDVI per county:
var means = ndvi_maine.reduceRegions({
  collection: maine,
  reducer: ee.Reducer.mean(),
  scale: 30,
});

// selecting only desired columns and transforming it into a feature object:
means = ee.Feature(means.select(['CntyFips', 'mean']));

means = means.setMulti({date: '2018-01-01'});   // trying to add properties

print(means);


texts: converting a list of shapely geometry to numpy array


texts: Are there any portable versions of QGIS?

texts:
I am looking for a portable version of QGIS.
Do you know if there are any portable version of QGIS?


texts: How to make sure that vector and raster data align correctly in ArcGIS?

texts:
I have a shapefile of European states (open-access from here: http://tapiquen-sig.jimdo.com) and DEM data downloaded from Copernicus.

Europe states vector is in GCS_WGS_1984
(Datum: D_WGS_1984,

if projected to ETRS89_ETRS_LAEA datum: D_ETRS_1989)

Europe DEM tiles (raster) are in ETRS89_ETRS_LAEA (Datum: D_ETRS_1989)

If I overlay both files in ArcGIS, from large zoom they look ok:

(ArcGIS on-the-fly reprojection)

but when I zoom in, I see there is a shift of ~200-400 m (I have not noticed any trend going for example south to north, i.e. increasing or decreasing discrepancy).

Is it my raster, or vector data which is not correctly aligned?

Datums in reprojected data looks the same: D_ETRS_1989

I have plotted the basemap layer, and it seems that my raster (DEM, ETRS89_ETRS_LAEA) is better aligned with satelite imagery basemap.

Here is output using Effect tool, DEM is black. When transparent, it is obvious that aligns to basemap better that vector shp (red).

The reprojection of Europe shp from WGS84 to ETRS89 did not help, neither to open reprojected vector in a new dataframe.

How do I "objectively" assure that vector and raster data will align in all the time? (without for example manual shift?)

My fear is, that I can spot that raster and vector data do not align due to shore and islands edge effect. But what if I will have the same shift in my WGS84 GPS point data, and I will just overlook it because I have now way to verify where the data come from?


texts: Is it possible to add an ArcMap map to a report in Access?

texts:
I am using Access to create printable maps where I can change parameters (border colors) on the fly...

Currently, I am exporting the maps to JPEG then pulling them in.  I would like to import MXD files so I can easily modify either the maps or frames on the fly.


texts: How to extract information from a tiff through the use of a shapefile?

texts:
I am having trouble in using the data from the South America amazon from a shapefile. What I was asked to do was to take the data from the loss of forest in a format where it is possible to make sense of it. But I am unable to pass the information to tif format, that I would use in Rstudio to make analysis with the temperature data that I have.

By my understanding a tif mask with the Amazonian deforestation is what i need, but I am open to other suggestions.

link to shapefile: http://www.imazongeo.org.br/downloads/imazon_sad_desmatamento_2017-01_amazonia.zip

ps: I have tried using QGIS to convert, but the dots are lost in the conversion, and I am unable to do it in R* too.


texts: Cartodb map style changes when clearing view

texts:
I'm losing the custom settings of my map when I perform simple filtering analysis and "Create Dataset from Query".

Here's what I'm doing:

I've created a map where I've specified my color scheme, updated the Info Window with the fields I want visible and update the legend. I then Publish the map.
In the DataView, I filter my data and "Create Dataset from Query"
I return to the MapView and the map retains my customization of the Info Windo, color scheme and legend.
There is an option to "Create Dataset from Query" or to "Clear View" in the MapView.  If "Clear View", I lose all my custom settings (Color scheme, no fields selected in Info Window and legend disappears)

I thought when I customized my map then Published it, every time I would "Clear View" the map would return to what I've saved?


texts: Appending Token to ArcGIS API for JavaScript Geocoder Widget

texts:
I am building a web map application using the ArcGIS API for JavaScript.  The application is using the GeoCoder widget from ESRI, but I need to link it to an outside, secured resource.  I am able to use the outside server to create a token from my credentials, but I don't know how to append the token information to the GeoCoder widget.

I am looking for the correct way to add the token information into my geocoding requests.  I have been scouring the forums here and on GeoNet, and the information I find points mostly to using either the IdentifyManager or using a proxy.  The outside server lets me create a token that is valid for 1 year, so I would prefer to do that right now and set up the IdentifyManager or proxy later.

My site is HTML/JavaScript based on an IIS (6.1) server, and I was not planning on making it an ASP site.

I would prefer if I just take my credentials and create a token to quickly insert into my application rather than taking the time to set up the other methods.

Is there a way to take the web-generated token and pass it into the GeoCoder widget?


texts:
I'm creating a map in ArcMap 10.6 and I'm having trouble exporting it to a PDF.

Everything exports fine to PDF except for some SPOT imagery that I've loaded through a WMS Server. The imagery displays in the data frame in ArcMap but, when I export to PDF, the imagery does not get exported although all other layers export perfectly.

I've also loaded the mxd into ArcMap 10.3 and the same problem occurs.

Has anyone found a way to export WMS layers to pdf through ArcMap 10.6 or 10.3?


texts:
I recently upgraded my Geoserver from 2.1.1 to 2.17.1
Since I have done that, one of my layer show unwanted lines in his polygones :

This is what is show by a GetMap WMS request on the layer
The layer is from a PostGIS table filled with a shapefile read by FME.
In the previous Geoserver, the layer is OK.
In QGIS, shapefile is OK and data from table is OK too.
I tried to republish the layer and recreate the table without success


texts: is there a way to manually stretch a polygon?

texts:
I have this polygon that I want to manually fit and resize(without entering affine numbers, just resizing it with my mouse).
Is there a way to do it?


texts: Interpolation in QGIS shows wrong output result

texts:
I have done interpolation (IDW and Nearest Neighbour) in QGIS. But interestingly the output shows the wrong result. The area where I have a huge number of points results in a lower Z value (it is supposed to have a higher Z value).
Is it because of the overlapping points?

Another doubt is that kriging algorithm does not show the predicted value (i.e. predicted value is not in the Z value range -160 to 80), though the output curve is meaningful. The prediction shows values in the range -120 to -70 or so. What might be the reason?


texts: Import Foreign Schema Fails with "schema postgis does not exist"

texts:
I am connecting to a Postgres 9.4 foreign server instance where Postgis 2.1 was installed using CREATE EXTENSION postgis WITH SCHEMA postgis from a Postgres 9.6 Beta where Postgis 2.3dev was installed with CREATE EXTENSION postgis in the public schema.

When running:

IMPORT FOREIGN SCHEMA data_schema FROM SERVER f_server INTO data_schema;

I receive the following error:

[2016-08-07 11:23:36] [3F000] ERROR: schema "postgis" does not exist
  Where: importing foreign table "data_table"
[2016-08-07 11:23:36] [25P02]

WHERE data_table is defined as:

CREATE TABLE data_schema.data_table(
  id bigserial PRIMARY KEY,
  geom geometry(MultiPolygon,4326)
);
CREATE INDEX ON data_schema.data_table USING GIST (geom);

I am aware that some postgis functions have schema qualification issues, as demonstrated by this postgis ticket; however that doesn't seem to be the issue I am encountering since I merely have a geometry column and a GIST index.

Tried importing the postgis schema from the foreign server as well, but that errors out with:

IMPORT FOREIGN SCHEMA postgis FROM SERVER f_server INTO postgis
[2016-08-07 11:43:08] [42704] ERROR: type "postgis.geometry" does not exist
  Where: importing foreign table "raster_columns"

Foreign Server Definition:

CREATE SERVER f_server
FOREIGN DATA WRAPPER postgres_fdw
OPTIONS (host 'localhost', port '5433', dbname 'f_db');

CREATE USER MAPPING FOR postgres
SERVER f_server
OPTIONS (user 'postgres', password 'password');


texts: message: "string.substitute could not find key \"loadingState\" in template" esri javascript 3.23

texts:
I am trying to add a measure widget from the esri samples and I get the following error.

I am using boostrap.css and I am wondering if this has an effect on widget loading.

I am also wondering if these tags conflict with each other. I have trouble getting widgets to show in my map div.

I toggle the classes to change the page width and show a side bar menu.

I saw a similar question with no answer.


texts: Applying merge policies when merging features in same layer using ArcGIS Desktop?

texts:
I am trying to figure out how to merge line features within the same feature class such that particular fields will be summed. For example, in the following picture I have two lines that I wish to sum the "TRVLTIM" or travel time field, such that the merged line has a value of 0.055737.

I have tried setting up a domain with the "Merge Policy" set to "Sum Values" and applying it to this field, but as described on the following ArcGIS help page, the "Merge" button in the editor toolbar ignores merge policies.

Is it possible to prevent the editor merge tool from ignoring domains?

I want to repeat this process for several roads. I am aware that geoprocessing tools exist that can perform this sum, but it would be impractical to recreate this massive roadway feature class every time I want to merge.


texts: Iterator (Iterate Feature Selection) running loop in undesired path using ArcMap ModelBuilder?

texts:
I am building a model using Iterate Feature Selection. I want to run the Iterator -> select by location -> Calculate field, and after finishing this loop (#3 in the image below), I will continue to the next process (#4 in the image below).

But when I run the model, it go to the next process until the end, then come to start the next loop, repeatedly.


texts: Calculate volume of space defined by multiple overlapping 3D polygons

texts:
I have thousands of 3D polygons scattered in space but they mostly overlap with many of their neighbors in big extent.

Is there any way how to compute volume they define effectively?

I would love to have precise volume, but some relatively good approximation will do too.

My solutions so far:

v.patch() (GRASS) - so far, I was able to just patch faces and kernels, but it does not form closed 3D polygon after looong computation.
ST_3DUnion() (PostGIS) - unfortunately, the function is not an aggregate function (it accepts only two simple arguments). I have wrote function calling ST_3DUnion iteratively, unioning polygons one by one. But obviously, computation is slow and basically stops on polygon 100.
Substracting polygons from master volume (PostGIS) - creating one big polygon from bbox of all my polygons and then using ST_3DDifference while iteratively substracting every single polygon. What I am left with is big polygon with a hole with volume of my interest. Never tried - probably same issue as 2).
Overlap with 3D point grid (PostGIS) - I could create 3D points not far away each other in such way they form regular grid which covers bbox of all polygons. Finally, I would create spatial index on this layer and count points falling inside any polygon (operator &&& and ST_3DIntesects() could help a lot). But it is extremely memory demanding and gives "only" estimation of the volume. Kinda rasterization by 3D raster.

Obviously 1), 2) and 3) are suffering from exponential increase of faces/points after processing several polygons. Maybe, there is a way which does not require so expensive operations? Especially when I am interested only in volume.

Just to explain, I am trying to implement relatively simple Local Convex Hull algorithm in 3D. I have solved 3D hulling in GRASS (see my previous question) since ST_ConvexHull() does not work in 3D.


texts: How can i programmaticly get the WKT from WKID?

texts: Is there some way (in C#/silverlight) to get the WKT from some web service with just the WKID?


texts: DEM data - Central Europe - noise in heights

texts:
I have downloaded DEM data from this page: http://www.viewfinderpanoramas.org/dem3.html. I have written simple program in C++ that takes "bricks" and build resulting height map.
Hovewer, I have problem with result. I am using Mercator projection and take every "gps pixel" from DEM data and recalculate its position to real pixels (same for border). It works fine (or at least it seams to).

As you can see on image, there is something weird. That "mess" in lower part is incorrect and I dont know, if data are incorrect, if I am using them incorrect (wrong bricks ?) or something in my program is wrong.

EDIT:

Here is algorithm for reading "bricks" heights:

//iterate all bricks and set its starting coordinates
double startLon = 19;
double startLat = 48;

///-------------------------
int tileSize = 1201;
double step = 1.0 / tileSize;

for (int y = 0; y < tileSize; y++)
{
   for (int x = 0; x < tileSize; x++)
   {

    short height = brick[x + (tileSize - 1 - y) * tileSize];
        //tileSize - 1 - y -> reverse brick coordinates

    height = swap_endian<short>(height);

    if (height == -32768)
    {
        height = 0;
    }

    lon = (startLon) + x * step;
    lat = (startLat - 1) + y * step;

    int pixelX = 0;
    int pixelY = 0;
    this->helper->ConvertCoordinates(lon, lat, &pixelX, &pixelY);

    if (pixelX < 0) {  continue;}
    if (pixelY < 0) {  continue;}
    if (pixelX >= this->width) {  continue;}
    if (pixelY >= this->height) {  continue;}

        //write it to map
        //if multiple height in same pixel, take only last one
    this->heightMap[pixelX + pixelY * this->width] = height;
    }
}


texts: How to split lines at every intersection from a single geometry using PostGIS?

texts:
The query I tried is this:
Create table test1 as
  Select
    row_number() over() as ID,
    (st_dump(st_split(input.geom, blade.geom))).geom as geom
  from roads input
  join roads blade on st_intersects(input.geom, blade.geom);

The error that pops up is like this:
ERROR:  Splitter line has linear intersection with input
SQL state: XX000


texts: Point style with Python MapScript does not work

texts:
Software details: CentOS 7, Python 3.6, MapServer v 7.2.1 is compiled from source.

I use Python MapScript to build WMS, it works well for styling Lines and Polygons, but not for Points.

For example, the Line code:

style = mapscript.styleObj(classObj)
style.width = 10
style.color = mapscript.colorObj(255, 0, 0)

or the Polygon code:

style = mapscript.styleObj(classObj)
style.outlinewidth = 10
style.outlinecolor = mapscript.colorObj(255, 0, 0)

draws features as expected, while the Point code:

style = mapscript.styleObj(classObj)
style.symbolname = 'circle'
style.size = 10
style.color = mapscript.colorObj(255, 0, 0)

draws nothing.

Here is output of layer.convertToString() for Point layer. It looks like valid mapfile section:

LAYER
    CONNECTION \"dbname='db' host='db' user='user' password='1'\"
    CONNECTIONTYPE POSTGIS
    DATA \"value FROM (SELECT name, value FROM points WHERE ST_Intersects(value, !BOX!)) AS subquery USING UNIQUE id USING srid = 3857\"
    LABELITEM \"name\"
    NAME \"points\"
    STATUS ON
    TYPE POINT
    UNITS METERS
    CLASS
    NAME \"points\"
    STYLE
        COLOR 255 0 0
        SIZE 10
        SYMBOL \"circle\"
    END # STYLE
    END # CLASS
END # LAYER

The 'circle' symbol is defined in symbolset. There are no errors in log files. I can't figure out what is wrong. Does anyone have such experience?


texts: mxd.SaveACopy fails to save file using Python script

texts:
I have been struggling to build a model using ModelBuilder and Python scripts. At the end of the model, I would like to write a script which saves the MXD file using the actual current time of the process.
Hence, I have used the following script, but it does not work, displaying "failed to save document" error.
# Import arcpy and datetime modules
import arcpy
import datetime

#Setting variables
mxd = arcpy.mapping.MapDocument("CURRENT")
now = datetime.datetime.now()
TIME = now.strftime("%H:%M:%S")
mxdName = "Analysis" + str(TIME)+ ".mxd"

#saving Script
mxd.saveACopy(%GDBpath% + mxdName)

So, I decided to simplify the script. The result is as followed:
# Import arcpy and datetime modules
import arcpy
import arcpy.mapping as map

#Setting Variables
mxdName = "Analysis" + ".mxd"
mxd = arcpy.mapping.MapDocument("CURRENT")

mxd.saveACopy(r"C:\Users\SOHCO\Desktop\New folder" + mxdName)

But I still receive "failed to save the document" error.
Could anyone help me out?
Note: The %GDPpath% (from the first script) refers to a parse path tool in order to apply the directory that the user defines for the tool. (the image is attached at the end of the question)

picture of the error:


texts: Loop to intersect points with polygons in R for each individual

texts:
I have a shapefile with points for 6 individuals and another shapefile with 6 polygons (one per individual). Every individual has the same id in both files.

I want to intersect both layers in order to obtain, for each individual, the points within the corresponding polygon.
I've tried with only one individual, using the function raster::intersect() and it works perfectly. However, I will have to repeat this process for other animals so I would like to have a loop that could do this automatically (but I am not very familiar with loops).

I've tried several options, and so far this was the only one that worked partially:

for (i in unique(SLO_repr$id)) {
         for (j in unique(HR95_all$id)) {
                 if ((!is.na(SLO_repr$id[i])) == (!is.na(HR95_all$id[j]))) {
                         a <- raster::intersect(SLO_repr, HR95_all)
                         writeOGR(a, dsn= getwd(), paste0("points_withinHR", "_", i), driver = "ESRI Shapefile",overwrite_layer=TRUE)}
         }
 }

The problem with this one is that, in each shapefile, all points within each polygon are included. So, the raster::intersect() function worked, but it seems I didn't select the i and j properly...?

To build this loop, I followed this question: Loop to check multiple polygons overlap in r

I was also wondering if I really need to have a nested loop? Because all id's are the same, so in the condition, when it gets on i for SLO_repr$id, it should select the same value for HR95_all$id?

So I've tried this loop:

 for (i in unique(SLO_repr$id)) {
                 if ((!is.na(SLO_repr$id[i])) == (!is.na(HR95_all$id[i]))) {
                         a <- raster::intersect(SLO_repr[[i]], HR95_all)
                         writeOGR(a, dsn= getwd(), paste0("points_withinHR", "_", i), driver = "ESRI Shapefile",overwrite_layer=TRUE)}
         }

But I get this error: "Error in as.vector(y) : no method for coercing this S4 class to a vector"

Any ideas about this?


texts: Refer to zoom level in TileMill expressions

texts:
Is it possible to do something like this:

#ways[highway="trunk"] {
  line-color:#c99;
  line-width:@zoom/10;
}

Obviously that example doesn't work. Is there another way to do arithmetic expressions involving the zoom level?


texts:
I have found Gevel, a tool enabling the visualization of the content of the GiST indices in Postgres. In particular, I am interested to use it to visualize the RTree spatial index in PostGIS. Does it work with Postgres 2.0? Is there a simple way (not compiling from source) to get it to run (e.g., with postgis.app on Mac OSX?) Apparently, Gevel and the rtree-gist modules should now be part of the default Postgres distro as per the readme, but I cannot find any way to enable these.

EDIT: Following on from comments, if I have to build from source, do I have to build while building Postgis, or can I add gevel later.


texts: extract scientific layers from MODIS HDF Dataeset using Python GDAL

texts:
I have been trying to extract scientific data layers from MODIS HDF files but without success. Before I have been doing this type of task using Arcpy, but I want the script that I am currently working to be based only in GDAL and not in Arcpy so I can easily share it.
I used GDAL, but it seems that GDAL version in windows does not have support for HDF files (I browsed some forums about it). I tried to open the HDF file using gdal.Open() but it does not return anything, and I am pretty sure that my directory path to the HDF file is correct.
Then I tried pyHDF. So this packages allowed me to extract the scientific layers as 2D array, but I have a problem how to save it as GeoTIFF file. I am aware that for me to save a 2D array to a GeoTIFF file I need some reference, such as projection and raster properties, just like the GetGeoTransform() module in GDAL. In this step, I am not sure how to extract that same kind of information using pyHDF.


texts: What is the meaning of "MULTIPOLYGON EMPTY" Message?

texts:
I am trying to run ST_Intersection function.

SELECT ST_AsText(ST_Intersection(ST_Point(42.729123082000058,-96.619302355999992), ST_Point(42.72899715300008,-96.619410439999967))) FROM clu_08aug;

I am getting Resultset with the "MULTIPOLYGON EMPTY".

Can someone help me what does it means?

Dharmesh


texts:
Is there, with QGIS 3.x, any way to restrict the data added to a specific field of a GeoPackage layer?

I want to restrict the data in a field of an attribute table with something like that:

2letters(a-z)-3letters(a-z)-2numbers(0-9)


texts:
I need to export a 3d model made from a DEM with a satellite photo draped over it. I'm using qGIS 2.8.2 and NVIZ, and am having a few problems.

1) How do I get a full colour image draped over a DEM in NVIZ? When I import the original (it's a jpeg but I can convert it) into GRASS GIS it separates it into RGBA channels, and I can only overlay one of those over the DEM.

Any suggestion on how to overlay the optical image?

2) I've got a working DEM with a single channel image overlaid, but i go to export using one of the three options from the 'save image as' menu (PPM image, TIFF image, or Max resolution PPM), nothing happens.

Does anyone know how to get these functions to work?

Or if they're working, where do they save the files/what do are they named by default?

I'm ultimately doing this so I can import a 3D scene into Blender; if anyone's done this in a more straightforward way I'd be happy to have your suggestions.


texts: Add and customise attribute table to print composer

texts:
I've finished a plan of farmland constituting of a polygon called 'parcels' which contains critical information for the final/printed plan. I could simply label the polygon as shown below, but there's a lot of clutter.

I've added a new field in the 'parcels' layer which associates a number to each feature, just like the ID field when creating a new shapefile. The idea for this is to label the polygon with a single number, then adding and tweaking the attribute table in print composer. Below is the look I'm trying to replicate, except I need a header such as "Number, Parcel ID, Ac, Ha"

How do I get a table as shown in the last image into the print composer? The atlas generator seems to be going in the right direction.


texts:
I want to insert an image in a popup when I click on a layer. I'm successfully using Leaflet-pip but when I want to display an image I've always have undefined. My url is stored in a postgres database.

myLayer.on('click', function(a){
        var match = leafletPip.pointInLayer(a.latlng, myLayer, false);
        if (match.length) {
            popUpText = '';
            for (var i = 0; i < match.length; i++) {
                popUpText += '
';
                popUpText += match[i].feature.properties.name + '

';
                popUpText += '
' + match[i].feature.properties.photo + '

';
            popUpText += '';
            }
        }
        var popup = L.popup().setLatLng(a.latlng).setContent(popUpText).openOn(map);
});

My field photo looks like this

''

How can I solve this ?


texts: Get version of installed QGIS from Windows command line

texts:
Tell me please, how can I get current version of installed QGIS from command line (batch-file). QGIS installed with Osgeo4w-installer.
Registry search did not return any result.


texts: Unsupported operand error writing to file (Python, arcgis)

texts:
In the code below, I get a "TypeError: unsupported operand type(s) for +: 'NoneType and 'str'" on the f.write line. If I don't concatenate, it runs fine, but all the layer names run together.

I understand what the error means, but not why I'm getting a NoneType. I've tried casting lyr.name to str, same result. All I'm trying to do is write out a text file of layer names. Any idea what's causing this and how to get around it? I've tried testing for type is None, isinstance(None), no luck.Thanks.

    f = open(r"C:\ListOfLayers.txt", "w")
    count = 0
    lyrs = arcpy.mapping.ListLayers(mxd, "", df)

    for lyr in lyrs:
        f.write(str(s)) + "\n"
        count += 1

    f.write("\n \n There are " + str(count) + " layers.")

    f.close()


texts: Projecting shp from GK3 to UTM32 with BeTA2007

texts:
I've experimented with projections of shapes, used the official provided grid beta2007 for germany and noticed differences between transforming by double exporting over custom CRS in QGIS and ogr2ogr of around 0.8-1.0 metres...
Perhaps it's happening because of false strings... I don't see it :)

Exporting
I've got the following custom CRS DHDN GK3 with the grid (EPSG:31467):

  +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +units=m +no_defs +nadgrids=X:\gridshift\BETA2007.gsb +wktext

Exported the shape with this CRS, loaded it into QGIS and exported it again to an custom CRS ETRS89 / UTM zone 32N with leading 32 (EPSG:4647), this is also the project CRS, on-the-fly projection is enabled:

  +proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=32500000 +y_0=0 +ellps=GRS80 +units=m +no_defs +towgs84=0,0,0,0,0,0,0

Doing the same thing shorter with ogr2ogr

  ogr2ogr -f "ERSI Shapefile" utm32.shp gk3.shp -s_srs "+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +units=m +no_defs +nadgrids=X:\gridshift\BETA2007.gsb +wktext" -t_srs "+proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=32500000 +y_0=0 +ellps=GRS80 +units=m +no_defs +towgs84=0,0,0,0,0,0,0"

Why does this happen?


texts: Adding database connection to ArcSDE database?

texts:
I am unable to connect to the SDE database we have in house on my new CPU and not quite sure what's going on. I don't have much experience with SDE but seems like the connection should be quite simple. I can connect to the SDE database on another computer using a username and password as it was already set up. But with this new computer, no dice. Not sure if there is something else that needs to be installed on my new machine to make the first connection or if there is something else I am missing. IT is trying to help but they are stumped as well. Here are stats for both CPU:

New CPU with no SDE connection:
Windows 10
64-bit
32 GB RAM
Intel Core i7-6700 CPU
Running ArcGIS for Desktop 10.2.2

Old CPU with successful connection:
Windows 8
64-bit
8 GB RAM
Intel Xeon CPU E5-1620 v2
Running ArcGIS for Desktop 10.2.2


texts: Selecting many maximum points within polygons

texts:
I've got a point layer with three values such as high, age etc. and polygon layer with many forest area. What I want to do is to select the 100 maximum (highest) points in every polygon.

I can't do it normally because I've got something like 10,000,000 points and over 10,000 polygons.

I'm using ArcGIS 10.2.2. for Desktop Standard, QGIS 2.8.6 and global mapper.
I'm looking for an answer not in every software. I need just a solution. It could be even in different software such as Access.
I already have in point layer column with polygons ID. But it is just too many of them.


texts: Creating NDVI time series graphs from buffered point (vegetation plot) locations in Google Earth Engine

texts:
I am aware there are a lot of similar questions to this on here like

Time series analysis of NDVI using Google Earth Engine
Creating NDVI time series from single polygon of shapefile using Google Earth Engine

but I couldn't quite find one that helped me.
I am looking at tree death from drought and I have on ground survey data as polygons(plots) and points.
I am trying to get an NDVI time series graph for each buffered point (user defined in map window as geometry2) using Landsat 8 32 day composites. Id like to see if NDVI can show the decline in tree health (and death) over time at these point locations (with 100m buffer).
Ideally I would like a mean value for each poly or point(buffered 100m) for every L8 32day NDVI composite from Dec2017-Dec2020(present).
I keep getting the error:
"image" is not defined in this scope.
in , line 122
But the code doesn't have a line 122.
// Load Landsat 8 Collection 1 Tier 1 32-Day NDVI Composite input imagery.
var landsat8NDVI32d = ee.ImageCollection('LANDSAT/LC08/C01/T1_32DAY_NDVI')
    .filterDate('2017-12-01', '2020-12-01')
    .select('B[1]');

// Define regions of interest as 100M buffers around points.
var geometry = ee.Geometry.Point(geometry2).buffer(100);

// Reduce the region.
var meanDictionary = image.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: region.geometry(),
  scale: 30,
  maxPixels: 1e9
});

// Create a graph of the time-series.
var graph = ui.Chart.image.seriesByRegion({
  imageCollection: landsat8NDVI32d,
  regions: geometry,
  reducer: ee.Reducer.mean()
})
print(graph)


texts: Custom function for mosaic seams in Python GDAL

texts:
I was wondering if there is a way to implement a custom function for creating seams between the rasters that create a mosaic. More specifically, if I create a mosaic with gdal_merge or gdalwarp, the images will be laid on top of each other one after the other without specifying where the seam between the images should be. I was wondering how (with which specific module in gdal) would it be possible to implement a function in Python that creates a seam like the one made with Area Voronoi Diagrams with  Overlap described here: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.725.3922&rep=rep1&type=pdf.
I'm very new to GIS, however I struggle to find anything of help on the internet and I wonder if it needs to be implemented completely from scratch... Can you point me in the right direction?


texts: How to source the boundaries(latitude/longitude) information for US states, cities, Zip level

texts:
I am trying to color code the states, cities, zips differently based on certain parameters. As the user zooms in they should be able to see the different colors dynamically painted on the maps. Was able to get this at state level as the boundaries required are available. However, struck at a point to get the boundaries for each city and zips of states. (http://www.usnaviguide.com/)

How to get these details? I have the enterprise licence of Google Maps as well. Do I need to procure any other license to get these boundaries so that I can paint them.

I had also gone through geo-coding service but I do not see any methods supporting to output the boundaries of area.

Thanks in Advance


texts: Localise qgis2web Leaflet output Measure tool

texts:
I just tested the measure tool in the leaflet output of the latest version of qgis2web plugin (2.28.0, QGIS 2.14.5LTR):

I've been trying to find a way to translate it to german and found that the /qgis2web/js/leaflet-measure.js contains a broad variety of languages, somewhere from line 6360... like

...
},{}],28:[function(require,module,exports){
// de.js
// German i18n translations

module.exports = {
  'measure': 'Messung',
  'measureDistancesAndAreas': 'Messung von Abständen und Flächen',
  'createNewMeasurement': 'Eine neue Messung durchführen',
  'startCreating': 'Führen Sie die Messung durch, indem Sie der Karte Punkte hinzufügen.',
  'finishMeasurement': 'Messung beenden',
  'lastPoint': 'Letzter Punkt',
...

way below (line 7048), I find some requirements, that might indicate that something is missing under /qgis2web/i18n/:

...
var i18n = new (require('i18n-2'))({
  devMode: false,
  locales: {
    'ca': require('./i18n/ca'),
    'cn': require('./i18n/cn'),
    'da': require('./i18n/da'),
    'de': require('./i18n/de'),
...

My last javascript experience is long ago, but are there some changes in the code that apply for a translation to german?


texts: Preventing rasterFromXYZ from changing raster extent?

texts:
I've created a new raster derived from a first one in which I set a specific extent.

e <- extent(-122, -18, -56, 14)
1st.raster <- crop(raw.raster, e)

This first raster is used as a mask for other analysis I've made. So the coords are the same through my modelling process.

My model result is a vector (Z values) with the same nrow of my 1st.raster. To produce a new raster with this vector I've extracted the coords from the 1st.raster.

coords <- na.omit(cbind(xyFromCell(current, 1:ncell(current)), values(current)))[,1:2]

2nd.raster <- rasterFromXYZ(cbind(coords, vec.Z))

Is expected that my e extent should be the same in both rasters, but that is not the case.

Here are some details of my real rasters.

> current
class       : RasterStack
dimensions  : **1680**, 2496, 4193280, 5  (nrow, ncol, ncell, nlayers)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -122, -18, -56, 14  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
names       : bio02, bio03, bio10, bio14, bio16
min values  :     9,    15,   -68,     0,     0
max values  :   211,    95,   303,   697,  3450

    > suit_c
class       : RasterBrick
dimensions  : **1680**, 1924, 3232320, 7  (nrow, ncol, ncell, nlayers)
resolution  : 0.04166667, 0.04166667  (x, y)
**extent      : -109.4583, -29.29167**, -56, 14  (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : in memory
names       :        SP1,        SP2,        SP3,        SP4,        SP5,        SP6,        SP7
min values  : 0.07558069, 0.07239924, 0.03424645, 0.07757321, 0.09786071, 0.09050093, 0.08286063
max values  :  0.8428620,  0.8791401,  0.8592408,  0.8562305,  0.8667494,  0.8393829,  0.8629798

I do need to create a new raster by combining the layers of this two rasters into one. As expected the different extents are preventing addLayer to work accordingly.

> FinalRaster <- addLayer(current, PA_c)
Error in compareRaster(c(x, rasters)) : different extent
>

How to create a new raster keeping the original extent?


texts: Converting dbf into raster image using Python as way to add table to layout?

texts:
ArcPy does not have a function to add a table to a layout.

I was thinking that maybe taking that table (e.g. a dbf), there may be some python script command to turn it into a raster (e.g. do something like a screen capture on the table).

You could then add the screen capture raster image to the layout using arcpy.

Anyone know of a function in python, or a few lines of python code to do this?


texts: Updating new field from existing field using ArcGIS Desktop?

texts:
I need to update my new field/column (Property Type) with either "detached", "semi", or "terrace" depending on an existing column (Property Count):

where the count = 1 I need to populate Property Type with "detached"
where the count = 2 I need to populate Property Type with "semi"
where the count >=3 I need to populate Property Type with "terrace"

How would I write this in the field calculator?


texts: GRASS GIS 7 in QGIS 2.6.1, Windows 7 64bit - Missing dependency

texts:
I am using QGIS Brighton 2.6.1. Today I have installed Grass GIS 7.0.0 RC2 to try its topology functions in QGIS.

I have activated the provider, set the paths, but everytime I try to use any Grass GIS 7 command, I get the "Missing dependency.This algorithm cannot be run :-(" error. Grass in standalone is working just fine. Windows 7 64bit user as it stands in the title.


texts:
I have an existing GeoServer instance running (version 2.1.3). This GeoServer has multiple Workspaces and Stores serving a number of layers as WFS to an OpenLayers map.

I want to update my Stores every quarter with the most recent layers. The process goes like this:

1) Delete the GeoServer Store (for example Store: "3rdQ2012") under Workspace "A" that contains the previous quarters shapefile.
2) Create a new GeoServer Store (for example Store: "4thQ2012") under Workspace "A" for the new layer.
3) Edit my OpenLayers map project to point to the new layer.

For some reason when I do this GeoServer is not serving the new layer. When I look at the GeoServer log file I see this message for the affected layer:

"2012-12-21 16:15:33,965 WARN [geotools.gml] - Could find type for nameOfLayer in the schema, generating type from feature."

When I go to the GeoServer "Layer Preview" page I can successfully see and select "OpenLayers" under the common formats column GeoServer successfully opens a map window and displays the layer. But, when I select GML under the common formats column GeoServer opens a new window and displays this :

" This XML file does not appear to have any style information associated with it. The document tree is shown below.
http://schemas.opengis.net/wfs/1.0.0/OGC-exception.xsd">
java.lang.NullPointerException null
"

It appears that GeoServer is not seeing my new layer or perhaps something is not updating correctly after updating my stores and layers.

I have tried restarting my GeoServer service as well as doing restarts on the server. I have also deleted the workspace, stores, and layers multiple times in an attempt to force GeoServer to update but all to no avail.

Any suggestions would be appreciated. Thanks!


texts:
I am trying to create a custom projection in GeoServer, but I see nothing (from the logs, neither in the SRS_list) coming out.
The projection I'm trying to create is a Mercator for Mercury.

I start from the following proj4 string:

+proj=merc +a=2439700 +b=2439700 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs

Then, in Python, I use pyproj to convert it to WKT:

>>> from pyproj.crs import CRS
>>> merc = CRS('+proj=merc +a=2439700 +b=2439700 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs')
>>> merc.to_wkt()
'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",2439700,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Reference meridian",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CONVERSION["unknown",METHOD["Popular Visualisation Pseudo Mercator",ID["EPSG",1024]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'

Which I then append to my GeoServer's /usr/local/tomcat/webapps/ROOT/data/user_projections/epsg.properties as:

90000000=PROJCRS["Mercator / Mercury", \
    BASEGEOGCRS["unknown", \
        DATUM["unknown", \
            ELLIPSOID["unknown",2439700,0, \
                LENGTHUNIT["metre",1, \
                    ID["EPSG",9001] \
                ] \
            ] \
        ], \
        PRIMEM["Reference meridian",0, \
            ANGLEUNIT["degree",0.0174532925199433, \
                ID["EPSG",9122] \
            ] \
        ] \
    ], \
    CONVERSION["unknown", \
        METHOD["Popular Visualisation Pseudo Mercator", \
            ID["EPSG",1024] \
        ], \
        PARAMETER["Latitude of natural origin",0, \
            ANGLEUNIT["degree",0.0174532925199433], \
            ID["EPSG",8801] \
        ], \
        PARAMETER["Longitude of natural origin",0, \
            ANGLEUNIT["degree",0.0174532925199433], \
            ID["EPSG",8802] \
        ], \
        PARAMETER["False easting",0, \
            LENGTHUNIT["metre",1], \
            ID["EPSG",8806] \
        ], \
        PARAMETER["False northing",0, \
            LENGTHUNIT["metre",1], \
            ID["EPSG",8807] \
        ] \
    ], \
    CS[Cartesian,2], \
    AXIS["(E)",east,ORDER[1], \
        LENGTHUNIT["metre",1, \
            ID["EPSG",9001] \
        ] \
    ], \
    AXIS["(N)",north,ORDER[2], \
        LENGTHUNIT["metre",1, \
            ID["EPSG",9001] \
        ] \
    ] \
]

I then restart GeoServer but nothing comes out in the logs (no error) neither in the Demos/SRS_list.
It is like the projection/file is not read.

Clearly, I'm missing something obvious here...Any thoughts?


texts: QGIS Raster to Polygons produce invalid geometries

texts: I trying to convert a raster layer into polygons with QGIS. At first, it seems there is no problem with doing that, but when I'm checking the output for invalid geometries, it always gives  me some invalid geometries (Ring Self-Intersections as shown in the picture). Does anybody know, how to avoid this problem or a workaround for it?


texts: How to visualize layers from ArcGIS\QGIS project in their order

texts: A lot of times I've seen visualization of layer in their order (like on the pic). Could anyone recommend a tool to make the same with my layers (could be QGIS or ArcGIS) or at least how this visualization is called?


texts:
I need to use the roads to compare GPS positions to verify that they are no errors in the data.

I downloaded the map.osm from a site that have information of City of Rome and had it transformed to Spatial Lines in R.

src <- osmsource_file(file = "rome.osm")
bb <- center_bbox(41.9, 12.3, 1000, 1000)
rome <- get_osm(bb, src)
ways <- find(rome, way(tags(k == "highway")))
ways <- find_down(rome, way(ways))
ways <- subset(rome, ids = ways)
hw_lines <- as_sp(ways, "lines")

So with the Spatial Lines Object I tried to use  distm function but It's accused that the spatial object can't be used (not is in a suitable format).

There is another way using OSM to get information of only roads that I can compare with GPS data in R?


texts: Explode MULTILINESTRING into individual segments in PostGIS 1.5

texts:
I have a MULTILINESTRING in PostGIS consisting of several LINESTRING's. Consider each LINESTRING to have a set of POINT's

  P1, P2, ... PN

forming the line segments

  L1, L2, ... LN-1

where LX is the line segment connecting PX to PX+1.

Using PostGIS 1.5 how can I explode the MULTILINESTRING into individual line segments of each LINESTRING?


texts: User interaction with Python tool midway through script running?

texts:
I need to design a Python script that will take a CSV input, search for a matching polygon in our SDE data for each line, and then extract certain information about data that intersects with the polygon.

This is straightforward, except that about 5% of the input values will not be exact matches (e.g. a parcel number will have hyphens in the CSV, but the SDE value has no hyphens). I've got a few tests in there to handle these cases, but it will inevitably result in some situations where the script returns multiple possible candidate polygons and the user needs to be able to indicate which is correct. For example, the number may pull up Polygon A and Polygon B as possibilities, and it's up to the GIS analyst to visually examine to decide which is in the correct region.

If I was running the script in IDLE, then I could simply use raw_input to get that user input and proceed with the rest of the steps. Is there an equivalent mid-script method in ArcMap that accepts user input? I am only aware of the initial tool parameters, which wouldn't be useful for this issue.


texts: ESRI Collector app- Offline record collection limit

texts:
I am wanting to find out what the maximum number of records people have collected while being offline on the ESRI Collector app. I am operating on an Ipad with a 64 GB capacity.

Has anyone load tested the Collector app offline or know what the limit could be?


texts: QGIS - use atlas attribute with 'Like' command

texts:
Does anyone know how to use the new get Attribute function in conjunction with the "Like" command?

I have a table with an Attribute "Sub_Region" and want to show the features that match the atlas feature of the same name...

The syntax below works for "="
"Sub_Region" = attribute($atlasfeature, 'Region_Name')

However, some of my precincts cover two regions, so the attribute might be 'Central, Western'. I can manually show only those of interest by selecting:
"Sub_Region" like '%Western%' .... however I want to inset the Attribute function so that I don't have to manually change the query as I export the atlas but can't work out the correct syntax...

Any suggestions?


texts: Set CRS to numpy array with no geographical data using rasterio

texts:
I have 1 dimension numpy array that has no geographical information. The numpy array was originally satellite image that I have manipulated in order to get this 1 dimension array.
I have the crs and transform of this array.
However, until now, in order to have geographic data on this array I had to write it as tiff and only then I could get back the geographical data and do things such as crop the raster and more.
I know there is option to reproject an array with no geographical data, as mentioned here . I have tried to do this with my data as following, based on the doc:
src_shape = (886, 1138)
rows, cols = (886, 1138)
d = 1.0/240 # decimal degrees per pixel
# The following is equivalent to
# A(d, 0, -cols*d/2, 0, -d, rows*d/2).
src_transform = A.translation(-cols*d/2, rows*d/2) * A.scale(d, -d)
src_crs = {'init': 'EPSG:4326'}
source = img_pred

dst_shape = (886, 1138)

dst_transform = rasterio.transform.from_bounds(*bbox, width= bbox_size[1], height= bbox_size[0])
# this gets the transform which is in this case : Affine(0.00011618510158013974, 0.0, -44.75374, 0.0, -7.050087873462236e-05, -6.80845)

dst_crs = rasterio.crs.CRS.from_dict(init='epsg:4326')
destination = np.zeros(dst_shape, np.uint8)

reproject(
    img_pred,
    destination,
    src_transform=src_transform,
    src_crs=src_crs,
    dst_transform=dst_transform,
    dst_crs=dst_crs,
    resampling=Resampling.nearest)

Then I get new array with values 0:

So my question is : how do I get my array with the values reprojected correctly?
I assume that the problem might be in the definition of the original CRS, (src_crs), as my numpy array has no crs at all, but not sure if this is the reason.


texts:
I'm trying to create a df through a GeoJSON url, merge it with a df created out of a SQL query, and create a shapefile/GeoJSON from the output. The output has everything put the .prj file (in the instance of a shapefile output). I've tried to set the CRS later in the code but it doesn't stay. I've discovered at some point in the merge, it drops the CRS.

import geopandas as gpd
import pyodbc
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

pd.set_option('display.max_rows', 500)
pd.set_option('display.max_columns', 500)
pd.set_option('display.width', 1000)

url = 'http is private just for this example.geojson'
df = gpd.read_file(url)
df = df.rename(columns={'MAINADDRESS': 'Address'})

cnxn = pyodbc.connect( xxxx)

sql = sql query
sql1 = 2nd sql query

df1 = pd.read_sql(sql, cnxn)
df1a = pd.read_sql(sql1, cnxn)

df4b = pd.merge(df1a, df1, on='strap')
df4c = df.merge(df4b, on='Address', how='left')
gdf = gpd.GeoDataFrame(df4c, geometry='geometry')
print(gdf.crs)
gdf.to_file('rentals.shp')
gdf.to_file('rentals.geojson', driver='GeoJSON')

The print(gdf.crs) shows 'None'.


texts: Features lost when extruding polygons between two TINs

texts:
I am trying to transform a collection of adjacent simple polygons to a collection of simple polyhedrons. My approach is to create two constant-elevation rasters representing the upper and lower boundaries of the polyhedrons, convert these two rasters to TINs, and then use the Extrude Between tool to generate the polyhedrons. However, I find that some of the features are lost (see gaps in the figure) and there appears to be some weirdness going on in the edge definitions of the resulting polyhedrons.

If I am not doing something wrong, is there an alternative approach to vertically extruding polygons to create polyhedrons?


texts: Extract angle of valley walls on both sides

texts: I have a DEM of a the bedrock surface beneath a glacier and what I want to do is subtract the mean angle of the valley walls on the left and right side of the glacier along this line. I have a line that divides the glacier into a left and right side (the brown dotted points). Any suggestions on how to do it efficiently?


texts: reading .tif and .tif.aux.xml as a whole using geotools

texts:
I have sample.tif, sample.tif.aux.xml, sample.qml files in a directory. The .qml and .xml files contain classification data about .tif file. For example my .xml file is something like this:

<PAMDataset>
   <PAMRasterBand band="1">
    <GDALRasterAttributeTable>
     <FieldDefn index="0">
       <Name>Class_ID</Name>
       <Type>0</Type>
       <Usage>0</Usage>
     </FieldDefn>
  <FieldDefn index="1">
    <Name>Class_info</Name>
    <Type>2</Type>
    <Usage>0</Usage>
  </FieldDefn>
  <Row index="0">
    <F>1</F>
    <F>Water</F>
  </Row>
  <Row index="1">
    <F>2</F>
    <F>Vegetation</F>
  </Row>
  <Row index="2">
    <F>3</F>
    <F>Built-up</F>
  </Row>
  <Row index="3">
    <F>4</F>
    <F>Bare_soil</F>
  </Row>
    </GDALRasterAttributeTable>
  </PAMRasterBand>
</PAMDataset>

I want to read these files as a whole. I tried this way:

File file = new File("sample.tif");
AbstractGridFormat format = GridFormatFinder.findFormat(file);
GridCoverage2DReader reader = format.getReader(file);
GridCoverage2D coverage = reader.read(null);
List<Category> categories = coverage.getSampleDimension(0).getCategories();
for(Category c : categories)
{
   System.out.println(c);
}

According to the javadoc for category (http://docs.geotools.org/latest/javadocs/org/geotools/coverage/Category.html), I expect to have 4 categories for this file: (Water, Vegetation, Build-up, Bare_soil) but I found only one category with name "No data".

How can I read these files to correctly initialize categories collection?


texts: Open Source parametrized/multi-dimensional caching tile server

texts:
For a project I am looking at setting up a caching tileserver using open source components. Point data is arriving from an external source and gets stored in postgis. This data has a datetime and clientid field in addition to the geometry data.

I would like users (using an open layers webapp) to select a timeframe (and clientid) for which to visualise these points. But since data traffic and repeated map requests will be high, I fear that just doing this with a WMS server will slow things down. A caching tileserver on top of my WMS could solve this. But what I don't know is if there are any tileservers out there that are aware of additional dimensions (parameters) beyond x/y/z.

I am currently looking at mapnik/mod_tile (skipping the WMS step), mapserver/mapcache and geoserver/geowebcache but am not restricted to these at all. Custom development on top of FOSS components is ok too (but I hate re-inventing the wheel).

Does anyone have any recommendations?


texts: Unable to create NDVI image using Python gdal package

texts:
I have Python 3 and all GDAL bindings installed on my Windows machine, so that when I run this:

import numpy as np
from sys import argv
from osgeo import gdal, gdalconst

I get no error messages. What I want to achieve now is to create NDVI image from two Landsat 8 images. To implement this, I'm using this nice looking script - https://gist.github.com/celoyd/7443646 . I really like how simple it is, however, when I run it like so:

C:\LC81750202017091LGN00>python ndvi.py LC81750202017091LGN00_B4.TIF LC81750202017091LGN00_B5.TIF output-ndvi.TIF
ERROR 1: Deleting output-ndvi.TIF failed:
Permission denied
ERROR 1: TIFFResetField:output-ndvi.TIF: Could not find tag 273.
Traceback (most recent call last):
  File "ndvi.py", line 31, in <module>
    ndvi = (n - r)/(n + r)
MemoryError

I tried to debug it and see that the first error message

  ERROR 1: Deleting output-ndvi.TIF failed:

comes from these lines of code (19-23):

output = geotiff.Create(
   argv[3],
   red.RasterXSize, red.RasterYSize,
   1,
   gdal.GDT_UInt16)

How can I fix it?


texts: Extracting values at specific latitude,longitude from MODIS swath data

texts:
I am trying to determine the amount of precipitable water vapor (PWV), ozone and aerosols as a function of time over specific spots on the Earth, namely our astronomical observatories. To do this, I've already got some Python code using modapsclient which will download the twice daily MODIS Aqua and Terra MYDATML2 and MODATML2 products that cover the specific latitude and longitude I am interested in.

What I'm not sure about is how to extract the specific quantities I want such as the time the MODIS data were taken and PWV for the particular latitude and longitude position of my observatory to make them into a time series of values. The MYDATML2 products seem to contain 2D latitude and longitude grids of Cell_Along_Swath_5km and Cell_Across_Swath_5km so I guess this makes it swath data as opposed to tile or grid data ? The quantities I want such as Precipitable_Water_Infrared_ClearSky also seem to be against the Cell_Along_Swath_5km and Cell_Across_Swath_5km butI'm not sure how to get that PWV value at the specific lat,long I am interested in. Help please ?


texts: Spatial subset using a rectangular bounding box - geometry not in Mercator?

texts:
I want to restrict my feature collection (in red, below) to the conterminous U.S. I found this solution to create a bounding box:

var continentalUS = ee.Geometry.Rectangle(-127.18, 19.39, -62.75, 45);

However I'm perplexed by this because displaying "continentalUS" using Map.addLayer() shows curved latitude lines that are more "curved" further north. The background Google Map, as well as everything else to my understanding, is projected to Mercator before it's displayed, and also before operations in the script are performed. I tried converting "continentalUS" to a feature, and still no luck with the projection matching when added to the map. This curved rectangle doesn't subset the target feature collection as I'd like. Here's what I get when I use .filterBounds() with the above rectangle coordinates:

I want the 45th parallel to match up with the map and my shapefile, which is automatically projected by Google Earth Engine. Why isn't a geometry created in Google Earth Engine not in Mercator, like everything else?


texts: How to determine the closest feature to my points?

texts:
I am doing a research project on transportation for motorcycles. For each of the trips, I have the origin point and the final destination point. I also have another land-use layer with my points of interests (E.g. hospitals, schools, residential etc.). For each of these origin/destination points, I want to know where do they come from and where do they go. What should I do to answer this question?

I do not have the license for network analysis.


texts: How can I create a map of towns that fullfill several conditions in QGIS?

texts:
I successfully created a map of Senegal with several indicators. Now I would like to do a heatmap with the goal to show villages that respect those criteria:

1) high population density;

2) village or town (no city neither hamlet);

3) less then 10km from an electricity grid;

4) low land elevation;

5) no protected area.

To be more clear, I am joining the map I have created with the important indicators:

LEGEND:

Population density is the yellow -> red raster;
Town are stars;
Electricity grid is the red line;
Land elevation is the white -> brown raster;
Protected areas are the green polygons.


texts:
I am using the Refactor fields in the Processing Toolbox to change order of attribute columns. 3 of the fields have red text in the Expression column and when I Run the changes I get an error:

Algorithm Refactor fields starting...

Parser error in expression "Site name": syntax error, unexpected COLUMN_REF, expecting $end

How do I fix this??


texts: Lon-Lat to Simple Features (sfg and sfc) in R

texts:
How can I convert lon-lat points to simple features (sfg), and then put them in a simple feature collection (sfc)?

Here is a MWE that doesn't work but is the closest I have gotten to.

library(data.table)
library(sf)
# The DT data.table is the data I have (but 10,000s of rows, each row is a point)
DT <- data.table(
    place=c("Finland", "Canada", "Tanzania", "Bolivia", "France"),
    longitude=c(27.472918, -90.476303, 34.679950, -65.691146, 4.533465),
    latitude=c(63.293001, 54.239631, -2.855123, -13.795272, 48.603949),
    crs="+proj=longlat +datum=WGS84")
DT[, rowid:=1:.N]
# The following two rows do not work
DT[, place.sfg:=st_point(x=c(longitude, latitude), dim="XY"), by=rowid]
places.sfc <- st_sfc(DT[, place.sfg], crs=DT[, crs])
# This should result in five points, which it doesn't
plot(places.sfc)

I'm trying to learn Simple Features (which is why I do not want to use library sp), and later need to run st_buffer on the sfc.

Maybe better to create the sfc directly, without an sfg per point?

I use data.table for speed-reasons (10,000s of thousands of points that are also analysed without geographical aspects).

I think I need an sfc of sfg-points, and not a MULTIPOINT-sfg.


texts: Python Error; 'Charmap codec can't encode character '\ufeff' in position 0'

texts: Please help with this error when trying to use MMQGIS plugin.


texts: ArcObjects 10/10.1 Get rid of Invalid layer symbol in TOC

texts: I've written an add-in (VS2010, C#) that finds specific invalid feature layers in ArcMap and fixes them.  After the layer is fixed, the red invalid symbol still appears in the TOC.  An activeview.refresh doesn't take care of it.  How does one programatically get that layer entry to appear as valid without manual intervention?


texts: Dividing points and polyline into segments in ArcGIS

texts:
I have a question about measuring distances for points in separate segments. As you can see in the attached picture, the green line is my polyline which I want to divide into four segments A, B,C,D. I want to measure the distance of each point in a particular segment to the polyline. I don't want the distance to overlap in the sense that a point in segment A should not take its measurement in segment B or C or D. A point in segment A, for example, should take its measurement at any place on the polyline that falls in segment A only.
I must add that the the points and the polyline are separate features. Also, the red lines are just to show you how I want my segments to be.
So first of all, I want to know how I can divide the polyline and the points into segments. Then I want to measure the distance of each point in a particular segment to the polyline.
Just to remind that, the polyline and the points are two separate features.
I am using ArcGIS Desktop (ArcMap 10.8)


texts: Using ArcMap tool to create a weighted points in each village

texts:
I have 3 layers:

layer of districts; (shapefile)
layer of cadaster(village) inside each district; (shapefile)
Geodatabase that containes info about each cadaster with latitude and longitute, number of houses, number of people inside each house and more info.

What I need to do is, for each cadastre I need to create a point that make the sum of all houses inside this cadastre, and I need to display on top of each point, the number of houses and the number of habitants.

I tried the aggregate tool, statistics tool, and the spatial analysis tools, but it didn't get me the right result that I really want.


texts: How do I check / set the ObjectID to user-controlled?

texts:
Where do I look to check how a feature class ObjectID is assigned? My understanding of the topic is there are two ways to populate the objectID field: let ArcSDE can supply it using a sequence or the user can supply an objectID.

I'm not looking to have a discussion about best practices and the pro's and con's of each approach.

I just need to know how I can check how the objectid for a feature class is populated and, depending on the what I find, set it to a user-defined objected.

This is ArcSDE 10.1 recently upgraded to 10.4


texts: I have a CSV file with the coordinates of the vertices of my polygons (in the picture they correspond to "left", "right", "top", "bottom"). Actually is a grid of polygons. How can I import the file in QGIS as polygons and not as points?


texts: QGIS print composer: can I select features whose attributes to show?

texts:
I am using the atlas to generate a number of maps.  Each map is of a particular zone represented by a multiline vector layer.  I have another layer with point data and each map has an attribute table with the point data.  It seems that I have just two choices:

I can have all the points in the table
or I can have just those points that are visible on the map

I want a third possibility:  just those points associated with the zone. i.e. an option to filter points relative to the 'current' feature -- something like "zone" = atlas.zone_name.


texts: Locating outliers in LAS Dataset with ArcGIS?

texts:
According to the statistics and symbology of my LASD I have outliers. See below:

How can I locate them in ArcGIS 10.x and do further analysis with them? I cannot query the dataset or select point? So how can I extract them?


texts: Avoid adding new generated intermediate layers into ArcMap during Python script running?

texts:
I am running my script to do geoprocessing in the python window of ArcMap.

Since the geoprocessing includes lots of intermediate files, and these files were added in ArcMap at the same time of the processing.

I thought it is a waste of time to display all these intermediate files in ArcMap.

Can anyone help to suggest a way to only operate on these intermediate files in background and without adding them to ArcMap?


texts: Flatten/smooth/reduce noise in SRTM elevation maps

texts:
Is it possible to "flatten" some areas in an elevation map using GRASS?
For example I have an elevation map of Iraq. I'd like to smooth out the mostly flat areas while still keeping things like rivers (tigris & euphrates), hills and mountains.


texts: Relative hyperlink paths in ArcMap 10 Map Document Properties?

texts:
Is there a way to use relative paths (instead of absolute paths) in the "Hyperlink base" field in Map Document Properties for an ArcMap 10 map document?

(I have "Store relative pathnames to data sources" checked, but somehow I still can't use relative pathnames in Hyperlink base)


texts: How to increase the number of the adjacent cells of a focal cell?

texts:
I want to recover the adjacent cells around a focal cell, I could do it by using the function adjacent form the package « raster »,
cells_adj <-  adjacent(Habitat[[1]], cells=cells_xy, directions=neigh[nei], pairs=FALSE)
where neigh <- c(4,8, 16) this function gives 4,8 or 16  cells through the argument directions, I want to include more cells (24, 48, 80, 120, and 168). Is there any option to do that?


texts: How to get image coordinate and corresponding location coordinate of a geo-referenced image in andriod

texts: In my android app, I have a geo-referenced image which I am going to use as map. So I need to divide it into tiles. For that I want its top left and bottom right image coordinate(x,y) with location coordinate of corresponding image coordinates using java code which is supported in android. Can anyone help me with the solution.