texts: Summing multiple years of rain data using Google Earth Engine?

texts:
I am trying to determine the what areas tend to get rain in a specific interval of time over multiple years. What I am doing is taking the filter date (in this case June 15th through July 15th) and outputting that data into its own layer. Where I am running into trouble is when I have multiple years of data, summing that data gives me a completely inaccurate output layer. What i'm trying to do is "stack" the years on top of one another on the map. See below for examples.

var dataset = ee.ImageCollection('OREGONSTATE/PRISM/AN81d')
var precip = dataset.select('ppt');
var precipitationVis = {
  min: 0.0,
  max: 50.0,
  palette: ['1621A2', 'white', 'cyan', 'green', 'yellow', 'orange', 'red'],
};
Map.setCenter(-86.55, 47.71, 7);

var precip2012 = precip.filterDate("2014-06-15", "2014-07-15");
var precip2013 = precip.filterDate("2015-06-15", "2015-07-15");

var total2012 = precip2012.reduce(ee.Reducer.sum());
var total2013 = precip2013.reduce(ee.Reducer.sum());

var final = total2012.add(total2013);
Map.addLayer(final, precipitationVis,
    'Final Layer');

Map.addLayer(precip.filterDate('2012-06-25', '2012-07-15'), precipitationVis, '2012');
Map.addLayer(precip.filterDate('2013-06-25', '2013-07-15'), precipitationVis, '2013');

This is the output for 2012

This is the output for 2013

And this is the output when I use the Reducer.sum function

So as you can see, the 'layered' final image is complete trash. I'm not sure where the numbers are coming from or how it is getting wild outputs (the minimum output value for the final layer is "100" when there at most in the previous layers were 4-5 mm of rain.)

What I want to be able to do is put 20 years of precip data for that june-july range onto one layer where it accumulates in each year if there was precipitation in that pixel. Does that make sense what i'm trying to do?

This is what the image looks like when I enable both the 2012 and 2013 layer at the same time:

And that's what I want the final product to look like, however when I do this just using the layers, the opacity goes to full after 3-4 layers.

So how would I compile each layer (2012 & 2013) into a final layer?


texts:
I sometimes have to write out the length along a line for every couple of hundred meters or kilometers. This applies to pipes, when the construction team wants to know where along the pipe they are or talk about.
So far I have done this in a very unefficient way:

dissolving line segments (if there are multiple along the way, and there usually are)
split them one piece at the time at specific distances using Editor --> Split
Convert endpoints to point features
Generate list of numbers in Excel to paste in a field in the point FC
label the points without symbol (or with a cross symbol as below)

This works, but there has to be more suitable ways of doing this. I'm not specifically looking for a programmatic/Python solution (even though it is welcome), anything that is a tad less complicated that the outlined workflow will do.. :)


texts: Snapping line endpoints to a point feature class by adding a vertex

texts:
I have a line feature class whose endpoints should snap to the closest point in a points feature class.

The problem is, I do not want the lines to move, but rather a vertex must be added to the current endpoint and the line extended so that its new endpoint is now coincident with the point.

If I use the topology rule on the line feature class that Endpoint must be covered by the point feature class, the topology correction offered is to create a new point on the endpoint of the line, which I do not want. Using the Snap tool in the Editing Toolbox doesn't move anything, and anyway I do not want the positions of my lines to be changed.


texts: Publishing AutoCad *.dwg files into HTML?

texts:
I'm new using Autocad. What I currently have is a dwg file generated from AutoCad 2015, now, I want to publish this file somehow into an html (web-mapping?) page so I can grab it from there and embed it into my own PHP application.

I've seen that usually the dwg files are published to Autodesk Infrastructure Map Server, but it's expensive software and I want to know if there's any other solution for this, if possible, open source and free.


texts: Exporting high-resolution world map image from QGIS?

texts:
I'm trying to get some (very) large world map images (think 10k x 10k pixels etc), but can't seem to do it with QGIS.

Manually adding smaller images to a larger one works (albeit slowly) for continental landmasses, but not for various mid-ocean islands etc, as I've no reference points handy.

The "Save as image" command only saves what's on-screen at the moment, the "Print Composer" system ditto.

Ideally, I'd hit a button & QGIS would export the entire project to a single giant .tiff image, or perhaps a series of tiles.


texts: Updating/replacing raster values in one raster values from another in QGIS?

texts:
I have two rasters: one that is a land cover map (cropland data layer) and the other is a rasterized polygon representing cranberry fields that aren't represented in the land cover raster that I have.  I would like to replace the values of the land cover map with the value from the cranberry raster.  The values of the land cover raster under the cranberry field raster are variable (from among 8 categories).

I assume the raster calculator is used for this, but I have little experience with that tool.

Image of the cranberry field raster:

Image of the land cover raster with cranberry raster overlaying:


texts: Linking two different Data Driven Pages in different data frames

texts:
I have a layer of data driven pages that go through this entire area, but the area covered is so large that it would be impossible to make out the page numbers in the inset map in the legend. What I am hoping to do would be to have the inset map cycle between the N, E, W, and S extents and then have a dynamic text element just saying which quadrant the map extents is currently located within. Each of these data driven pages have an attribute value denoting which quadrant they fall in (hence the symbology).
Is there a way for the inset map to grab a value from the main map's data driven page and then show the corresponding quadrant or is there another way to do this?
I am using ArcMap 10.8.1.


texts: Non-infrastructure distributed vector processing in R

texts:
Is there a non infrastructure R option similar to disk.frame for geospatial vector distributed processing? I have come across the geospark package but am not sure if it requires spark to be installed.
Disk.frame is a local non-infrastructure install that allows for distributed processing similar to that of Dask, Spark, etc. Below is the description:
https://diskframe.com/

How do I manipulate tabular data that doesn’t fit into Random Access
Memory (RAM)?
Use {disk.frame}!
In a nutshell, {disk.frame} makes use of two simple ideas

split up a larger-than-RAM dataset into chunks and store each chunk
in a separate file inside a folder
provide a convenient API to
manipulate these chunks

{disk.frame} performs a similar role to distributed systems such as
Apache Spark, Python’s Dask, and Julia’s JuliaDB.jl for medium data
which are datasets that are too large for RAM but not quite large
enough to qualify as big data.

As far as use cases I am thinking of point-to-point, point-in-polygon, etc. type of analysis, specifically for R.


texts:
I'm totally new to the field of GIS.

I want to create some GIS-models and publish them on my website. Therefore the models should be connected to my database to constantly have the most recent data.

At the moment I'm working with QGIS Desktop (3.4.2) and created some models, based on some PostGIS data served in a remote database and I also managed to export those models as a leaflet-based website. But it's only a static website with the data  requested from the database at that moment.

Now I tried to export / create a dynamic model, that I can publish and that holds a connection to the database, so that I can constantly request the most recent data from the database. But I could not find a method to achieve this.

I installed GeoServer, QGIS server with a QWC2 and also had a look at exporting a QGIS model as a standalone Python-script, but couldn't find a solution. I thought it might be possible to use QGIS as a kind of GIS-IDE (comparable with R / Shiny if you are familiar with that).

Is it possible at all to create such a dynamic GIS model for the web?

Is it possible with QGIS or do I need another tool to accomplish this?

Edit:

A quick example might be the best explanation: At the moment I can visualize all people going from A to B. To narrow the results down to a certain date or date range i need to set appropriate where-clauses in the SQL query. Once having a result I could export his as a website. But this would be static for the set date (range). To get results for another date or date range I would have to export the site once again for these new constraints. And this I would like to automate so that I can just pass arguments somehow to get the results for the set parameters.

Edit2:

After a lot of research i figured out that maybe Geoserver is the tool i need. I have the possibility to create SQL views, also based on parameters and Geoserver serves me the corresponding WMS-response. So far so good. But i still could not figure out how to export QGIS layers to Geoserver. There is this plugin called "Geoserver Explorer" for QGIS but as far as i found it it's only available for QGIS<3.x but i am working with QGIS 3.4.
Is there any alternative to "Geoserver Explorer" for QGIS >= 3.0 or how do you publish QGIS projects on Geoserver?

Edit3:
I also managed to publish a QGIS project via QGIS Server and it might be visible through QGIS Web Client 2 but then i ran into a deadend regarding passing parameters to the project hosted on QGIS server.


texts: I am trying to open *.hdf file from ASTER in QGIS and it can not recognize CRS. Does anybody knows solution to this?


texts: GeoServer layers containing cities in each country, given that all the cities reside in the same table

texts:
Using GeoServer, I need to provide a graphical representation of the distribution of some metrics on a geographical basis. The geographical regions will be coloured and labelled according to the value of that metric.

I would like to have a "world" layer (containing all the continents and only continents), a "Europe" layer (containing only countries from Europe) and a "France" layer (containing all the cities in France). So to speak, "a tree of layers containing the immediate children", given that these entities (continents/countries/cities) reside in three distinct tables with one-to-many relationships defined.

I haven't found any "filtering" option when defining a layer from a PostGIS store, so a "cities" layer can only be based on the whole cities table (not on the cities in country X). Is there something I missed ?

P.S.
I would rather have distinct layers than having one big layer and filtering with every request, because I think the caching feature can only refer to layers. Is this assumption correct ?


texts: zonal statistics sum on raster with same values

texts:
I was wondering if anyone knows what exactly zonal statistics sum does if all my input value raster have the same values.

For e.g. Lets say I have a finer resolution raster with 10 pixel each having the value of 10. I also have another coarser raster so that it perfectly overlaps the finer raster. I use the coarser raster as a zone to calculate the sum of all the 10 grids using zonal statistic sum. The resultant raster should have a value of 100 if I run the zonal statistics sum. But instead it has a value of 10.

I do not understand. Is it something that I am doing wrong or is it that unique values are required to do a zonal statistics sum?


texts: Using intersect in QGIS Python gives error "no geometry"

texts:
I am using the QGIS Python Console to automate some bulk processing. I want to create a Voronoi polygon layer from point data (solved this part) and then from a folder of shapefiles clip the Voronoi polygon layer with any of the shapefiles that intersect it. but I am getting the error:
AttributeError: 'str' object has no attribute 'geometry'

I am used to ArcPy and have really been struggling with understanding QGIS's Python. I feel like the cookbook and API resources are a little overwhelming for me (a beginner python-er)

Any advice?

import os
from qgis.gui import *
from qgis.core import *
import processing

#------------------------------------------------------
# User Variables:

#path to data location
soil_path = "D:\\soil\\soildata"
shp_path = "D:\\soil\\field_boundaries"
print ("Soil data: " + soil_path)
print ("Field Boundaries: " + shp_path)

#-------------------------------------------------------

def get_data(path, ending):
    data_list = []

    for root, dirs, files in os.walk(path):

        for data in sorted(files):
            if data.endswith(ending):
                data_list.append(os.path.join(root, data))

    return data_list
    print (str(len(data_list)) + " shapefiles found")

def clip(vor_poly):
    clipped_data = []
    shp_data = get_data(shp_path, ".shp")
    print (str(len(shp_data)) + " field boundaries found")
    for f in shp_data:
        print (f)
        print (vor_poly)
        if f.geometry().intersects(vor_poly.geometry()):

            OUTPUT = os.path.splitext(polygons)[0] + os.path.basename(f) + ".shp"

            processing.runalg('qgis:clip', vor_poly, f, OUTPUT)

            print ("clipped " + OUTPUT)

            clipped_data.append(OUTPUT)
    return clipped_data

def qgisprocess(soilpts):

    #voronoi polygons with 0 size buffer

    print("Voronoi Polygons")

    for f in soilpts:
        split_path = os.path.splitext(f)[0]
        inlayer = f
        outlayer = split_path + "_vor_poly.shp"

        processing.run("qgis:voronoipolygons", { 'INPUT':inlayer, 'BUFFER': 0, 'OUTPUT': outlayer})

        clippedpolys = clip(outlayer)

def main():
    #list of data
    soil_data = get_data(soil_path, ".shp")

    print("running QGIS processes")
    qgisprocess(soil_data)

if __name__ == '__console__':
    main()


texts: OpenLayers 3 - Feature Popup without Jquery/Bootstrap

texts: How do you create tooltips (context) on OpenLayers 3 map, using standard Javascript and flat CSS.  (No extra libraries)


texts: Dynamically set field information in popup template discussion (Automate the process)

texts:
I want to create a certain function to streamline some of my customization work.

I work with a web app that's built using ArcGIS JavaScript API which contains a few layers.

Each layer is configured to an identify task that uses a popup template to display attribute information from the layer. Now, there are lots of attributes PER layer that I manually configure up in the popup template declaration.

Ex.

var popupContent = {
title: "Business ID#: {UNIQUEID_1}",
fieldInfos: [
{fieldName: "BUSNAME_1", visible: true, label:"Business Name: "},
{fieldName: "Standard_Name", visible: true, label:"Standard Name: "},
{fieldName: "LATITUDE_1", visible: true, label:"Latitude: "},
{fieldName: "LONGITUDE_1", visible: true, label:"Longitude: "},
{fieldName: "INFO_EMP_1", visible: true, label:"Employment Info: "},
{fieldName: "Est_Employment", visible: true, label:"Estimated Employment: "},
{fieldName: "INFO_SAL_1", visible: true, label:"Salary Info: "},
{fieldName: "Est_Revenue", visible: true, label:"Estimated Revenue: "},
{fieldName: "PRMSIC4_1", visible: true, label:"Primary SIC: "},
{fieldName: "NAICS_1", visible: true, label:"NAICS: "},
{fieldName: "NAICS4_1", visible: true, label:"NAICS 4 Code: "},
{fieldName: "NAICS4DESC_1", visible: true, label:"NAICS 4 Desc.: "},
{fieldName: "NAICS2_1", visible: true, label:"NAICS 2 Code: "},
{fieldName: "NAICS2DESC_1", visible: true, label:"NAICS 2 Desc.: "},
{fieldName: "Use_Type", visible: true, label:"Use Type: "},
{fieldName: "UseSubcategory", visible: true, label:"Use Subcategory: "},
{fieldName: "UseCategory", visible: true, label:"Use Category: "}
]
};

This is one of the shorter templates that i use, but as you can see, it gets pretty tedious to set these manually. Once in a while is okay, but the layers are constantly getting fields added/removed.

I was wondering if any of you had any clever ideas as to how I could automate this.

One possible way I was thinking is to have a table in my ArcGIS server that acts as a lookup table for the fields.
Table structure:

The column field names in this table will be the "Layername_fieldname" and "LayerName_label" of the layers in my web app (from the layers' corresponding attribute tables).

Ex.

|| Layer1_fieldname || Layer1_label || Layer2_fieldname || Layer2_label ||
--------------------------------------------------------------------------
|| BusinessLoc      || Biz Location:|| CountyBoundary   || County Name: ||
||                  ||              ||                  ||              ||

Rows can be added or deleted as needed to each column as my web app layers grow or shrink in attribute size.

Basically, I want to know if this is a viable option or if there are other methods out there to perform this task. If this is a good solution, is it possible to query the field lookup table and use those values in my popup template constructor.


texts: Extract complete pixel values inside a geometry

texts:
I'm trying to write an earth engine script which extracts the pixel band values inside a polygon. (geometry)  Only for pixels completely within the polygon.  The polygons will be quite small. (field size)

I have looked through the docs, and searched for this, but am unable to find a solution.  Have only found statistical reductions.

// get the dates
var startDate = ee.Date.fromYMD(2018,1,1);
var endDate = ee.Date.fromYMD(2018,2,1);

// Create a geometry representing an export region.
var myGeometry = ee.Geometry.Polygon(
[[[-61.15985870361328, -32.94904673848672],
      [-61.152305603027344, -32.944797148619905],
      [-61.15659713745117, -32.9395388987784],
      [-61.16432189941406, -32.94378874137744]]]);

var landsatCollection = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR");

// Filter image collection to dateRange of interest
var landsatDateCollection = ee.ImageCollection(landsatCollection
.filterDate(startDate,endDate));

//Filter to AOI.
var landsatAOI = ee.ImageCollection(landsatDateCollection
.filterBounds(myGeometry));

// Filter to image with least clouds
var landsatBestImage = ee.Image(landsatAOI
.sort('CLOUD_COVER')
.first());

// Load a landsat image and select three bands.
var landsatDownload = ee.Image(landsatBestImage)
.select(['B4', 'B5']);

/////////
//Added code that was improved after Jon's comments
/////////

// generate a new image containing lat/lon of the pixel and reproject it to Landsat projection
var coordsImage = ee.Image.pixelLonLat().reproject(landsatDownload.projection());

var joinedImage = coordsImage.addBands(landsatDownload);

var valuesList = joinedImage.reduceRegion({
  reducer: ee.Reducer.toList(4),
  geometry: myGeometry
}).values().get(0);

print(valuesList);

Jon's inspiration on the MODIS script allowed the modification above.  The valuesList now contains all the pixel information.

To see how to Export the data in valuesList to a CSV see

Earth Engine convert list with coordinates and values into a feature collection for export


texts: Is the ground track of any GNSS satellite constant? (or does it drift?)

texts:
The from Wikipedia Global Positioning System pageone can see that the GPS satellites have repeating ground track

Orbiting at an altitude of approximately 20,200 km (12,600 mi);
orbital radius of approximately 26,600 km (16,500 mi),[86] each SV
makes two complete orbits each sidereal day, repeating the same ground
track each day

Is the ground track of the satellites belonging to the other GNSS systems (Galileo, Glonass, Beidou) always the same?
Because of nonidealities, such as solar radiation pressure, the satellite orbits, and thus, the ground tracks have errors. Will the ground track be corrected back to the "average ground track" somehow, or is it possible that the ground track drifts with time? (during months or years)


texts:
I can't really find the words to explain what I'd like to achieve so I attached a picture.

The top picture was made of two layers, the layer with circles was made/styled using the instructions from an answer to a question that I posted a few days ago (How to scale with radius attribute). The other layer was created by converting that same layer using the Points to Paths plugin. The bottom image is what I want to achieve. For background info, the data is about wind information from typhoon Haiyan that struck the Philippines last November 2013. I want to get/generate the storm swath (the lower picture from the image above) that kind of look like the screenshot that I took from USAID (image below).

How can I achieve this or is this possible? Thanks in advance.


texts: How to change a numpy array's dimensions

texts:
I'm using GDAL Python API to read a raster into a NumPy array, it will return a array's shape like [bands, rows, cols], if we want to use OpencCV to deal with this array, it will cause some problem, because OpenCV read an image into array' shape like [rows, cols, channels], how can I transfer an array[bands, rows, cols] into an array[rows, cols, channels] and transfer it as fast as I can, I have been using a way to do this job, but it cost a little much time, especially when the data is huge. Downbelow is my way:

import gdal
ds = gdal.open('some raster')
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount
ds_array = ds.ReadAsArray()
gdal2cv2 = np.zeros((rows, cols, bands))
for i in range(len(ds_array)):
    matrix = ds_array[i, :, :]
    gdal2cv2[:, :, i] = matrix


texts:
I am having trouble with using Python together with my gp.Calculatefield tool in my python script, and I am not sure why..
Heres the code (UPDATED)

import arcgisscripting, os, sys
myTempWorkspace = r"D:\ScheduledTasks\Temp\TempWorkspace_HyMS_Licences.gdb"
#workspace set up
gp = arcgisscripting.create()
gp.toolbox = "management"
gp.Workspace  = myTempWorkspace
gp.overwriteoutput = 1

codeblock = """\
def calField(Location):
    fldLocation = !LICENCE_2_LICENCE_0_LICENCE_LOCATION!

    if fldLocation <> " ":
        return !LICENCE_2_ROAD_NAME! + " (" + !LICENCE_2_LICENCE_0_LICENCE_LOCATION! + ")"
    else:
        return !LICENCE_2_ROAD_NAME!
"""

try:
    #gp.CalculateField_management(myTempWorkspace + os.sep + "LICENCE_3", "LOCATION", "[DWWEB_ADDRESS1] & \" (\" & [DWWEB_LOCATION] &\")\"","VB","")
    gp.CalculateField_management(myTempWorkspace + os.sep + "LICENCE_3", "LOCATION", "calField(Location)","PYTHON","codeblock")

    print "Calculate Field - Success"
except:
    print "Calculate Field - Failure"

I am getting the following error in the pythonWin IDE

>>> Unhandled exception while debugging...
Traceback (most recent call last):
  File "D:\temp\CalcFieldTest.py", line 21, in ?
    gp.CalculateField_management(myTempWorkspace + os.sep + "LICENCE_3", "LOCATION", "Location(newField)","PYTHON","codeblock")
RuntimeError:
exceptions.NameError: name 'Location' is not defined
Failed to execute (CalculateField).
>>> Calculate Field - Failure

Now I am not sure why I am getting the error, I think I have done everything I should within the codeblock and on the gp.calculatefield line.

I'm using ArcGIS 9.3.1
So any help is appreciated.


texts: Failing to build PostGIS on Mac OS X 10.6.8 (Snow Leopard)

texts:
I'm getting the same problem expressed here

Failing to build PostGIS on OS X 10.7 Lion

But in Mac OS X Snow Leopard. Things I have checked include:

PostgreSQL compiled to several architectures. I used the binary from http://www.postgresql.org/download/macosx/, and tested the dynlib objects with the tool file
PROJ and GEOS compiled for i386 and x86_64. I used the frameworks from http://www.kyngchaos.com
GDAL compiled from sources following KyngChaos instructions. I compiled it instead of using the provided framework because I'm a GDAL developer, and I'm going to work with it.
Trying to compile PostGIS, because I'm a PostGIS developer too, and I don't want to use the binaries. It has been configured correctly. I used this: https://gist.github.com/2002495

But I'm still getting the error: ld: symbol(s) not found for architecture i386

Other things I've tried:

Compile ONLY for arch i386. Not working. PostgreSQL expects PostGIS compiled for both architectures. I get the error ld: symbol(s) not found for architecture x86_64.
Compile ONLY for arch x86_64. Not working, for the same reason. I get the error ld: symbol(s) not found for architecture i386

I could provide configure and makefile logs, if needed, but the problem is this. It's a problem known in PostGIS lists, but I still haven't found a solution...

Thanks in advance, and best regards


texts: Batch Convert MapInfo .TAB Files to .KML

texts:
Possible Duplicate:
  ogr2ogr convert all shapefiles in a directory

As the title suggests, i have a few hundred .TAB files to convert to .KML.  This is likely to be a fairly irregular thing so FME (while brilliant!) is not really an option, and FW Tools doesn't appear to offer batch conversion, without possibly scripting the solution.  There there isnt a way in MapInfo/MapBasic to do this either.

Does anyone have any ideas?

Thanks in advance

AM


texts: Assigning custom symbol/overlay in categories of unique values?

texts:
Where I work at we update certain layers of data every month. I've recently been given access to upload data to Pictometry but I'm having an issue with some symbology that hasn't been explained very well.

In Pictometry explorer, the Zoning layer that I'm trying to update displays with an overlay of PUD over areas that are assigned one zoning purpose but also have PUD zoning:

In the CONNECTExplorer display you can see that they're all together in one symbology and the PUD shows up as an overlay. I thought that maybe there were two separate layers uploaded to go into one but you cannot do that with Pictometry.

The fields for this shapefile look like this:

My boss mentioned using the Overlay field in a query of some sort to pull them out into another layer but I already have completed that task.

And in ArcMap I have this in my table of contents:

Notice that I do not have a PUD field in ArcMap, this is the issue that I'm getting to.

I'm trying to figure out how to get PUD to be in my legend and to display as I want without having to use the option of Unique Values To Many Fields. I know that you can just make a layer package or something but that won't upload to the Pictometry software as it only accepts shapefiles.
I'm truly stumped on this and would love to just get this figured out for future reference.


texts: Perform Multiply (Times) function on several reclassified rasters using Python

texts: I have several reclassified rasters (0 = snow cover, 1 = no snow) and I would like to use the multiply function on all rasters that have the same spatial extent.  In ArcGIS 10 this function can be performed on only two rasters. Is it possible to do it on multiple rasters in Python?


texts: Displaying multiple sea level scenarios?

texts:
I'm trying to display a number of sea level scenarios to discuss landscape differences across time. I've created a CTM of the region now I need to display sea level at a number of reduced levels (e.g. -50m).

I've already tried extract by attribute using "value <=-50" which produces a number of error problems which seem to be related to the GDB, but I can't get round them. Could extract by mask be better?

Does anyone have a simpler or better way? In an ideal scenario I would like to display the final model with terrain changes displayed for the terrestrial data but the sea can be extracted as no data and displayed blue. However it can simply be a blue green display just showing coastline if that is the best way using ArcMap 10.3, Windows 8.


texts:
I am using QGIS 2.18.7 on my laptop.  I can add a new shapelayer, but when I stop editing, the shapes that I have entered disappear.  I have tried point, line and polygon shapes, they all appear to allow creation of points/line/polygons, but when I stop editing, it all disappears.

On my desktop machine it works correctly.

Is this a configuration or installation issue?


texts: Automating execution of Python console in QGIS?

texts:
I want to automate the execution of Python console in QGIS. When QGIS start, Python console starts  at the same time, and Python predetermined code wants to be executed.

Currently, in order to perform this behavior, the following procedures are required.

Plugins -> Python Console
Show editor -> Open file -> Run script

By setting up shortcut ("Setting -> Configure shortcuts"), starting of Python console (1) is possible.

However, subsequent processing (2) will be manual operation.

Is there any method of automating python console execution through processing (1) and (2)?

Python version: 2.7.4 (default, Apr  6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)]
QGIS version: 2.0.1-Dufour


texts: Just wondering if anyone knows when and how to do the transition to GDA2020 in QGIS? I assume it will be as simple as selecting the CRS but it doesn't appear to be available as a CRS yet. Currently using QGIS 2.14.3 - Essen.


texts:
I have read the "technical" documentation for the TIGER geographic datasets and they fail to describe the way the geometry is encoded. The document claims that the shape files encode the geometry of linear features, like roads, but the documentation shows no explanation of this. For example, the beginning of the spec for the all lines files is as follows:

As you see, it is just a bunch of meta data, not geometry. Not only that, the fields themselves are not explained. For example, MTFCC is described as being a "feature class", but the possible values or feature classes are not described! At the end of the file there is the following:

The last two fields appear to be some kind of key possibly linking to some other file (???) which has nodes in it (???) I have no idea because the spec does not further describe these two fields.

I have two questions:

(A) Where can I find an actual REAL spec that actually describes each field in the TIGER shape files?

(B) Where is the geometry for linear features located? (As in what exact file and field and what is the spec for those fields)?


texts: OD matrix "Error: For forward traversal, the AtJunction must be the FromEdge's ToJunction"

texts:
My OD Matrix creates this "Error: For forward traversal, the AtJunction must be the FromEdge's ToJunction".

Any suggestionshow I can fix this problem.


texts: Pixel Type issue for raster related tools in ArcGIS 10.5

texts:
I have found that the pixel type resulted from running a raster related tool on ArcToolBox on ArcGIS 10.5 was different from that from running exactly the same tool on ArcGIS 10.1 and ArcGIS 10.2.

For example, the "Feature to Raster" conversion tool on ArcGIS 10.5 generates a raster with pixel type of unsigned integer with a pixel depth of 1 Bit while the same tool on ArcGIS 10.1 achieves a raster output with pixel type of unsigned integer with a pixel depth of 8 Bit.

How does that happen given with the same input feature class? Is there a way to configure the pixel type and depth as a default setting on ArcGIS 10.5?


texts: Extent of shp files with R for overlapping

texts:
I have 2 shapefiles: data (spatial point dataframe) and polys (polygons dataframe). I want to do an overlap but is seems that it does not work.

Here are data and polys:

> data
class       : SpatialPointsDataFrame
features    : 12527
extent      : 10.20075, 20.6108, 54.08669, 57.75905  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
variables   : 3
names       :    timestamp_pretty,     timestamp,     imo
min values  : 01/04/2006 00:00:55, 1143849655232, 9048392
max values  : 30/04/2006 23:59:36, 1146441576823, 9191541
> polys
class       : SpatialPolygonsDataFrame
features    : 436375
extent      : 4210000, 5441000, 3395000, 4813000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
variables   : 2
names       :    Id, Count
min values  :     0,     0
max values  : 99999,     9

to manage the overlapping, I use

proj4string(data) <- proj4string(polys)                            # to confirm the same reference system
inside <- !is.na(over(data, as(polys, "SpatialPolygons")))    # overlapping shape file and data

and then mean(inside) to check the average of points in polys.

But there is nothing happening, mean is always 0. I used this many times before and it always work, I guess it does not work not because the extent of the 2 sph files are  different.

Is there a way to edit this?


texts: How to Bind attribute fields to combobox?

texts:
I'd like to list attribute fields in a combobox. there are two combobox . cmblayers and cmbfields . cmblayers list map layers . user select a layer in cmblayer and attribute fields of selected layer list in cmbfield. i wrote cmblayer but i can't write cmbfield code.How can i write cmbfield using xaml or code behind?
 cmblayer code :

<ComboBox x:Name="cmblayers"  ItemsSource ="{Binding ElementName=Map, Path=Layers}"    >

                     <ComboBox.ItemTemplate>
                                    <DataTemplate>
                                        <TextBlock Text="{Binding LayerInfo.Name}"/>
                                    </DataTemplate>
                                </ComboBox.ItemTemplate>

</ComboBox>

i found the below code . but it just bind to specify layer. i couldn't bind cmblayers to cmbfields

<ComboBox x:Name="cmbFields" ItemsSource="{Binding ElementName=Map, Path=Layers[IncidentsLayer].LayerInfo.Fields}">
    <ComboBox.ItemTemplate>
        <DataTemplate>
            <TextBlock Text="{Binding Name}"/>
        </DataTemplate>
    </ComboBox.ItemTemplate>
</ComboBox>


texts:
I want to create a WPS-chain in Geoserver.

Here's the HTTP-Post statement:

<?xml version="1.0" encoding="UTF-8"?>
<wps:Execute version="1.0.0" service="WPS" ...>
  <ows:Identifier>gs:Flst</ows:Identifier>
  <wps:DataInputs>
    <wps:Input>
      <ows:Identifier>FlstID</ows:Identifier>

       <wps:Reference mimeType="text/xml" xlink:href="http://geoserver/wps" method="POST">
        <wps:Body>
            <wps:Execute version="1.0.0" service="WPS">
            <ows:Identifier>gs:createnumber</ows:Identifier>
            <wps:DataInputs>
                <wps:Input>
                <ows:Identifier>plht</ows:Identifier>
                    <wps:Data>
                    <wps:LiteralData>testid</wps:LiteralData>
                    </wps:Data>
                </wps:Input>
            </wps:DataInputs>
            <wps:ResponseForm>
            <wps:RawDataOutput>
                <ows:Identifier>number</ows:Identifier>
            </wps:RawDataOutput>
            </wps:ResponseForm>
            </wps:Execute>
        </wps:Body>
       </wps:Reference>

     </wps:Input>
    <wps:Input>
      <ows:Identifier>FlstFeatures</ows:Identifier>
      <wps:Reference mimeType="text/xml" xlink:href="http://geoserver/wfs" method="POST">
        <wps:Body>
          <wfs:GetFeature service="WFS" version="1.0.0" outputFormat="GML2">
          <wfs:Query typeName="Test:nyc_census_blocks">
            <ogc:Filter>
            <ogc:PropertyIsEqualTo>
                <ogc:PropertyName>blkid</ogc:PropertyName>
                <Literal>360050028001007</Literal>
            </ogc:PropertyIsEqualTo>
            </ogc:Filter>
             </wfs:Query>
          </wfs:GetFeature>
        </wps:Body>
      </wps:Reference>
    </wps:Input>
  </wps:DataInputs>
  <wps:ResponseForm>
    <wps:RawDataOutput>
      <ows:Identifier>result</ows:Identifier>
    </wps:RawDataOutput>
  </wps:ResponseForm>
</wps:Execute>

The first WPS (gs:flst) has two Inputs:
1. A number, which is generated by a WPS called gs:createnumber
2. A SimpleFeatureCollection, which is retrieved from a wfs

The number from the first Input is the ID of the Feature in the WFS. How can I insert this number in the Filter-block?

I tried to use

<Filter>
<ogc:PropertyIsEqualTo>
            <ogc:PropertyName>blkid</ogc:PropertyName>
            <ows:Identifier>number</ows:Identifier>
</ogc:PropertyIsEqualTo>
</ogc:Filter>

but that didn't work (Exception org.geotools.filter.AttributeExpressionImpl cannot be cast to org.opengis.filter.Filter)

Thank's for helping!


texts:
I would like to read a raster layer from GRASS 7.8 into R 3.6.1 (on Ubuntu 18.04). Whenever I use the reasRAST method I get a cryptic message that simply says "no stars import yet" [sic]. The console capture below exemplifies:

> library(rgrass7)
Loading required package: XML
GRASS GIS interface loaded with GRASS version: (GRASS not running)
> library(magrittr)
> library(stars)
Loading required package: abind
Loading required package: sf
Linking to GEOS 3.7.1, GDAL 2.4.2, PROJ 5.2.0
>
> initGRASS(
+   gisBase = "/usr/lib/grass78" ,
+   home = tempdir(),
+   gisDbase = "/home/duque004/GRASSDATA/",
+   location = "global_igh_local",
+   mapset = "duque004",
+   override=TRUE)
gisdbase    /home/duque004/GRASSDATA/
location    global_igh_local
mapset      duque004
rows        61319
columns     160300
north       8600750
south       -6729000
west        -20037500
east        20037500
nsres       250
ewres       250
projection  NA
>
> r <- readRAST(c("bio_zones_02"), mapset = "BIOLOGY")
Error in .read_rast_non_plugin(vname = vname, NODATA = NODATA, driverFileExt = driverFileExt,  :
  no stars import yet
>

What could be causing this? I have used similar code countless times with other combinations of GRASS and R, but never saw something like this.


texts: Calculating average slope from line using QGIS?

texts:
I have a raster dem from which I've created a slope. Then I'd like to calculate the average and maximum slope on each of some vector lines. I have a vague memory that I have seen a plugin for this but right now I can't seem to find it?

Basically I'm looking for something like the Zonal Statistic Plugin but instead of a polygon, I want to use lines.


texts: Determining which layer in raster stack is first to meet condition using R?

texts:
I want to analyze a raster stack and determine the first layer in the stack that meets a condition for every pixel.

Here's a reproducible code example generating rasters and making a stack.

# Generate stack of random values
r1 <- r2 <- r3 <- r4 <-r5 <- r6 <- r7 <- raster(ncol=10, nrow=10)
r1[] <- rpois(ncell(r1), 1)
r2[] <- rpois(ncell(r2), 1)
r3[] <- rpois(ncell(r3), 1)
r4[] <- rpois(ncell(r4), 1)
r5[] <- rpois(ncell(r5), 1)
r6[] <- rpois(ncell(r6), 1)
r7[] <- rpois(ncell(r7), 1)

# stack them
s <- stack(r1,r2,r3,r4,r5,r6,r7)

I want to analyze this stack to determine which layer is the first to meet a condition like: "Value = 1 for 3 consecutive layers" or "2 of last 3 layers have value = 1", for example.

Can you think of any function to do this?

If I run rasterToPoints I can convert the stack into a matrix with columns representing each layer, this may make things easier.

Once I determine the correct layer corresponding to each pixel, I then want an output raster with that layer number as the value for each pixel.

A somewhat similar question has been previously posted as Counting layer number of raster stack with value under some conditions and write it into output layer?, though that's counting the number of layers that meet a condition.


texts:
I have a list of paths of all the shapefiles (around 50) to be merged. All of them are at the same projection and none of them is sized more than 110MB (shp+dbf). The combined size of merge (shp+dbf) is around 1.5GB. The arcpy code used is

arcpy.Merge_management(list_of_paths, output)

Problem:
Merging all these files somehow creates an empty file (after around 15 minutes) with an error message:

arcgisscripting.ExecuteError: ERROR 999998: Unexpected Error.
Failed to execute (Merge).


texts: X/Y Coordinates to Longitude/Latitude, package for .net core app

texts:
I have a task to convert some incoming geo data to Longitude/Latitude. The data elements in the json are identified as xCoordinate and yCoordinate.

However, the data in my sample does not look like any geographic x/y coordinate data I've looked at online, so I have a couple questions that I could use expert help on, please forgive the lack of knowledge regarding this data, I've never worked with geo data prior.

An example from my sample dataset:

xCoordinate: 0981829
yCoordinate: 0197525

Couple questions:

Do these look like valid X/Y coordinates? Online examples seem to be in decimal format.
I need suggestions for a package to import into a .net core application which converts X/Y data to Longitude/Latitude.


texts: Data organization strategies for single user desktop QGIS

texts:
I have been teaching myself GIS using QGIS for several weeks and I am looking for advice on how to structure my data in a way that is extensible and won't require me to constantly reorganize things.  I have been working with the following types of data: topography (DEMS, Hill shades, contours), hydrography, roads, trails, political boundaries and some associated/joined data files.

Currently my data is in a mis-mash of file types and naming conventions based on the source I got the data from.

My questions:

should I use a database, and if so, what which one: SQLite, SpatiaLite or Postgres?  Or should I simply use the file system?
are there some rules of thumb for naming/structuring the data I have outlined above?

I realize this may not be the most crisp way of asking the question, but I am looking for a wholistic approach to organizing my data, since it is very painful to reorganize things across multiple projects on a regular basis.  Looking for any hints/guidelines you can provide.  It seems like many people must use the same general data types and have ways of organizing it in a way that is fairly standard?


texts:
I'm trying to build a dictionary from a table view in ArcMap using SearchCursor. Here is the snippet of code I'm having issues with:

valueDi = dict([(key, val) for (key, val) in
             arcpy.da.SearchCursor
             (joinFc, [joinIdFld, joinValFld])])

This is the error I'm geting:

Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 323, in RunScript
    debugger.run(codeObject, __main__.__dict__, start_stepping=0)
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\__init__.py", line 60, in run
    _GetCurrentDebugger().run(cmd, globals,locals, start_stepping)
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\debugger.py", line 654, in run
    exec cmd in globals, locals
  File "N:\Common\Script\Drivers\PrepCleanScapesUpdate.py", line 1, in <module>
    import arcpy
ValueError: too many values to unpack

I've checked the brackets and parentheses against a snippet of code used in an answer to Using UpdateCursor for joined field in ArcPy?


texts:
I have a lot of stacked points/coordinates. With the point displacement function inside QGIS I can render these stacked dots into a circle.

But I like to organize these points into a grid like this example:

How can I achieve this?


texts: Using ModelBuilder to calculate multiple polygon areas into one table?

texts:
I am attempting to calculate the areas of a large number of shapefiles (all polygons) in ArcGIS 10.2.

All of the separate shapefiles are stored in the same group layer, and each contains one polygon that needs to have its area calculated. The actual .shp/.dbf, etc. files are all stored in the same folder.  It would be best if the areas could be tabulated into one .csv after being calculated. I have considered using the "Calculate Area" tool in the Spatial Statistics toolbox along with an iterator, although I am not proficient enough with model builder and iterators to get it to work. I can go through and do each shapefile's polygon individually, but there has to be a faster way.

Any hints?


texts: How to import HDF (L1G) files into ArcGIS for Desktop?

texts: I have downloaded satellite imagery from GLCF website, it has downloaded with .L1G file extension, I am not able to open this file in ArcGIS as well in Erdas, so please let me know how to open this.


texts: Raster-extract function - Area-Weighted values

texts:
I am working in R and want to calculate the value of a polygon. The value should consider the weights on each intersecting cell. When I try to run the "extract" function with a sample raster and polygon I get different weights that the ones I calculate manually, resulting in different final value.

Here is my sample code:

r <- raster(nrow=2, ncol=2, xmn=-180, xmx=60, ymn=-30, ymx=90)
r[] <- c(1,2,4,5)
s <- raster(xmn=-120, xmx=-40, ymn=20, ymx=60, nrow=1, ncol=1)
s.pl <- as(s, 'SpatialPolygons')
w <- raster::extract(r, s.pl, method="simple",weights=T, normalizeWeights=F)
mean.value <- raster::extract(r, s.pl, method="simple",weights=T, fun=mean)

The value I get is 2.14 but according to the actual weights of the cells it should be 2. Can it be because of the projection that is in lat/lon? But even when I assign projection in meters I have the same result. How can I get the value of 2 that I am interested in? I also tried with "resample" function but I get different results as well.

My final target is to create a new raster with different resolution and extents from the original one and assign the values based on the weights of the cells of the original raster that intersect with the cells of the new raster. But it seems that neither resample nor extract functions give the expected outcome.

Hope my question makes sense.


texts: Developing standalone gis applications using pyqgis

texts: I am trying to develop a standalone application using pyqgis but I have no clue on where to start and where to write the code to make it happen so if anyone can provide me a tutorial to getting started, it'll be a  huge benefit for me. Thanks.


texts: Merge different KMZ layers to a unique KMZ file, a bash script

texts:
I have different .kmz files, each one contains a .png and a .kml, let's say:
a.kmz, b.kmz, c.kmz.
I would like to use a bash script to obtain a unique .kmz file, let's say d.kmz, that contains all the previous layers.
I simply tried to compress the .kml files and call the result d.kmz, but it does not work.


texts: Changing units project wide from 'Millimeter' to 'Map Units'

texts:
All my units in my project are in millimeter. But I would like to switch  all units to map units at once without clicking every single signature. I also would like to have the possibility to directly switch between units for the whole project.
Is there a way?


texts:
I am a true novice when it comes to ArcGIS and computer programs like this one. Am taking a class and cannot get the Project function to work under the tool box.

The path I used was ArcTool box --> Data Management Tools --> Projections and Transformations--> Feature Path --> Project  --> (This gives me a blank box and when I double click on the OK button I get the following message:

Error 000735:  Input Dataset or Feature Class:  Value is required
Error 000735:  Output Dataset or Feature Class:  Value is required
Error 000735:  Output Coordinate System: Value is required

Can anyone tell me what I can do to remedy this?


texts: What's the accuracy of Web Mercator EPSG:3857?

texts:
Are there any references on the accuracy or (min/max/avg) error of Web Mercator when it comes to e.g. distance calculations?

For the sake of completeness, I haven't found any reference to the accuracy Web Mercator. In contrast, spatial computations in WGS-84 projection are, due to Charles Karney, highly accurate as described in http://arxiv.org/abs/1109.4448 This, however, requires Gnomonic projection of coordinates for point to line projections which is very costly in GIS software. An alternative is the use of de-facto standard Web Mercator which must be questioned if accuracy is sufficient.


texts:
I am looking for a way to clip image collection by a given polygon geometry with Python API,
The GEE JavaScript equivalent is simply ImageCollection.clip(geom)
(https://developers.google.com/earth-engine/image_visualization#clipping)
With the python API I get the error: 'ImageCollection' object has no attribute 'clip'.

I'm adding my code after using .mean() as a reducer
My code:
set the boundary polygon
my_boundary = ee.Geometry.Polygon([ [[-5, 40],[65, 40],[65, 60],[-5, 60],[-5, 60]]] ,proj=None)
#set the ImageCollection and Image
img_Collection = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_NO2').filterDate('2020-03-01', '2020-03-27').select('NO2_column_number_density')
my_Image=img_Collection.mean().clip(my_boundary)
#Create Ipyleaflet map
my_map = ipyleaflet.Map(zoom=7)
my_map.add_layer(layer=my_Image)
my_map
error:
AttributeError: 'Image' object has no attribute 'model_id'


texts:
I downloaded congressional district and county shapefiles from nhgis, i then ran an intersect geoprocess on the two which ended up deleting some of the congressional districts.

This geoprocess overwrote my original shapefiles so I deleted the files and re-downloaded the shapefiles from nghis, but when I opened them up arc the cd file was still missing the districts.

Not sure what is going on here, perhaps a problem with a lock file? I tried renaming the files and looking for a .lock file that I missed somewhere, but no luck. I also tried closing arc and restarting that to see if that would remove the lock but that didnt work either. I am not too familiar with .lock files and have never run into this problem before so I am pretty stumped on where to start looking for answers.

EDIT:
im now thinking it is a problem with the original data I am using, could someone confirm this? You can download the data from nhgis, search by geographic levels click all then click Congressional District (1987-1993, 100th-102nd Congress) (by State) and  then download it (you do need to have an account). Then see if you are seeing the same thing I am in arc and check the .dbf file to see if is missing arizona (statefp 04) congressional district 01, search 0401.


texts:
I use geoserver geowebcache for displaying tiles using google maps API. Data comes from Postgis. My problem is that I have to change geometry of some features (not very frequently).

I would like to achieve automation simillar to this from openstreetmap, where after source data change tiles that cover area of changed feature regenerate themselves automaticly. How to do this? Now I have to delete all tiles from geoserver and generate them once again. I do not use pure WMS as service is rather popular and WMS may work slowly...


texts:
I have recently started to learn how to use open source web maps, and have come across the following issue.

When using geoserver and geo explorer (part of the open-geo-suite), I seem to get an annoying white box around my geotiff (see pic) almost like its been re projected or something. The geodetics of the geotiff are WGS UTM30. Ive tried re projecting the image to different projections but im having zero luck..

Has anyone got any ideas?


texts: Looking to find (or build) a particular Python geoprocessing tool

texts:
There is a large area (a wildlife refuge) of nearly 200,00 acres that we have loads of spatial data for, from habitat management projects from the past several years. My boss (who isn't GIS savvy) would like to be able to choose a point or area within those 200,000 acres, and easily access all of the spatial data in it's vicinity, without having to navigate through an unwieldy number of layers.

I've put together a "master" mxd of all the layers, which looks like this:

What I'm looking for is a tool whereby my boss could input a latitude and longitude coordinate of interest, and it would return all the nearby spatial data (i.e. the points/lines/polygons it falls on or are within x distance).

Would it be possible to write a Python script that would do this?


texts: How to separate states in USA file using ArcGIS for Desktop?

texts: I have found a few different files containing the entire United States, with attributes including states. Is there a way to separate the states out into their own files (perhaps in one file geodatabase)? I later would like to use "Selection-> Select by State." While I can select each state individually and save them in one file geodatabase, I am wondering if there is a more efficient and quick way to do this. Also, any options I have found will not include both North Carolina and North Dakota, etc., because the state name is two words long. Instead I get one 'North.' I do not yet know Python.


texts: Create a new polygon from another one that is larger than the first by x meters

texts:
I have a set of polygons of different shapes and sizes in shapefile1.

What I need to do, is to create new polygons around the boundaries of the polygons in  shapefile1, larger by x meters.

Any suggestions?


texts: Create GPX track in R. writeOGR function

texts:
I have a gpx file where someone exported two different activities in two separate locations to a single track.  I need to split them and avoid getting the erroneous line in the track layer that connects the two places where they were using the GPS.  I haven't found a method to cut a track at a certain spot, so my approach was to read in the track_points layer, subset for the date I was interested in, and push those points to a spatiallinesdataframe, and write the new track.  I manage to write it without getting an error, but when I try to read it back again, it says there are no features in the track layer. I would be open to any solution that would split the track layer, and if anyone that is already answering this question knows how to write a 'tracks', 'track_points', and 'waypoints' layer to the same .gpx file, please include that in your answer as that will be my next step after solving this problem.

Here's what I did so far:

library(rgdal)
gpxfile<-"C:/Users/.../file.gpx"
trackpoints<-readOGR(gpxfile, layer="track_points")

library(lubridate)
trackpoints$time<-ymd_hms(trackpoints$time)
trackpoints<-trackpoints[trackpoints$time>'2016/03/16 12:00:59+00',] #subset for the date I'm interested in

coord<-as.data.frame(coordinates(trackpoints))
lin<-Line(coord)

lin1<-Lines(list(lin), ID="track")
lin2<-SpatialLines(list(lin1))

df<-data.frame(len=sapply(1:length(lin2), function(i) gLength(lin2[i,])))
rownames(df)<-sapply(1:length(lin2), function(i) lin2@lines[[i]]@ID)

lin3<-SpatialLinesDataFrame(lin2, data=df)

proj4string(lin3)<-CRS("+proj=longlat +datum=WGS84")
names(lin3)<-"name"

writeOGR(lin3,
     dsn="C:/Users/.../newtrack.gpx", layer="tracks", driver="GPX",
     dataset_options="GPX_USE_EXTENSIONS=yes", overwrite_layer = T)


texts: "Complex located in the street", due to messed up with CRS?

texts:
I need to combine two .shp files from two different sources.
First one is 20 Expanded parcels during 2012-2017 by MVP-CA, which I use as the base layer. The second one is a point layer of the Housing complexes of Beijing.
The problem is the CRS of two layers are different. EPSG:4326,WGS:84 for the first layer and

for the second layer. I have enabled the 'on the fly' option.
However, when I open two layers, all the complexes located in the street instead of parcels, where they are supposed to locate. I guess it's due to the contradiction of projection system, but I don't know how to fix it. Since I know directly change projection system can be problematic.

Update:
   I know there is an intentional deviation added to the GPS by the Chinese mapping authority to "keep secret of nation". Could this be the source of my problem?


texts:
I am trying to download time series data (2017-01-01 to 2018-12-31) of different optical indices (e.g. NDVI) from the COPERNICUS/S2_SR dataset. The chart I plot only ever contains few values per year and if I export the .csv I seem to have the complete data from Oct 2018 onward but most of 2018 and 2017 is missing.

I was working with the same code end of April/beginning of May and at that time 2018 data was more or less complete.

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


texts:
We have a folder accessible by QGIS server containing some projects to be served as WMS:

Projects are available in the qwc admin qui when creating a new theme:

There is a theme for each project:

I wonder why not all the themes are available in QWC2 theme section:

I tried deleting browser cache as well as ctrl-F5 (Firefox), but the themes are still missing.
Any advice one some kind of "refreshing" the QWC2?


texts: Postgis: ST_VoronoiPolygons gives TopologyException: Input geom 1 is invalid: Self-intersection

texts:
I have a problem that ST_VoronoiPolygons started to give me GEOSVoronoiDiagram: TopologyException: Input geom 1 is invalid: Self-intersection

The function that I use it for have worked for a long time for me.


texts: Finding index of max value in stack of rasters and reordering another raster stack based on index

texts:
I have two distinct tasks to complete- find the index of the maximum value in each cell of a stack of 12 rasters. There are 12 rasters because each raster is the mean monthly precipitation over 30 years. Then I want to reorder another raster stack of 360 layers, where each layer is the monthly precip over 30 years such that the resulting raster stack has cells that start at the index found to be the of the 12 mean monthly precip values and then continues with the rest of the monthly precip values.
So for example, I have monthly precip values Jan -Dec 1981, 1982 etc. However, I have found that the wettest month is July i.e. 7th month), so I want to reorder such for each pixel the stack starts from July 1981 until Dec 2020.
If I were to do this with a list of vector and not rasters, I would do
twleve_values<-c(11.630832  19.372524  16.788150   4.051670  14.857163 186.173774 420.738594 304.162227 161.079527  23.955036   8.945979   3.019485)
max_tweleve_value<-which.max(twelve_values) #this would result in 7 because the max value 420.738 is the 7th value

if(max_twelve_value > 1) other_360_values[c(max_twelve_value:(length(other_360_values) - (13-max_twelve_value)))] else other_360_values

I am stuck at the first task of finding the index of the wettest month in the stack of 12 rasters. I have tried
x<-stackApply(12_raster_stack, indices = rep(1,12), fun=which.max, na.rm = TRUE)

Error in FUN(newX[, i], ...) : unused argument (na.rm = TRUE) and
x<-calc(mean_monthly_ppt, fun= whiches.max)

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : cannot use this function
How do I move ahead?


texts: How to join repeated values to the same feature using QGIS?

texts:
I new to QGIS, and programming in general.

I'm trying to join a csv file with another file that is already in the project. However, the csv file has repeated values that should allocate to the same feature in the project, which doesn't happen. When I join the values only one of the repeated values appears; the others don't.

How do I solve this problem?


texts:
I'm a newbie and I need a little bit help with postgis query construction. My situation looks like this.
I have a table in postgresql with ways converted from bicycle routes (take from Open Street Map).
If there are 2 or more routes on the same way, the lines are overlapping on each other and have exactly the same geometry but have different tags (name, ref etc)
I would like to add another column with a counter for these duplicated lines. Let's say that I have 3 bicycle routes on the same way, this equals to 3 lines overlapping each other. These should get an extra tag with a number 1, 2 and 3. This should happend only if there are more than 2 overlapping lines. If there is only one, no number should be added.

I tried ST_Equals but can't get it to work....

Here's my code which is not exactly for this purpose but for merging lines with the same attributes. I would like to add the counters into it.

DROP TABLE  IF EXISTS  public.bicycle_merge;
CREATE TABLE public.bicycle_merge AS
SELECT ST_LineMerge(ST_Collect(way)) AS way,reversed,route, osmc_color,name, network, state, "instance:cycle", "relation:id", ref, description, distance
  FROM planet_osm_line
  WHERE route IN ('bicycle')
  GROUP BY   reversed, route, osmc_color,name, network, state, "instance:cycle", "relation:id", ref, description, distance


texts: How to close watermarks window in QGIS Openlayers plugin?

texts:
I'm using Quantum GIS Ver. 1.7.4. I have upgraded openlayers plugin 0.7 to 0.9. However when I open google satellite, a watermarks/copyright dialog box "Data Peta@2012 ..." appear in the middle of the screen. The close dialog box sign (X) appear on the top right of the dialog box but not active:

How can I close the dialog box?

Thanks.


texts: MapBasic Dialog Preserve command error: No active Dialog

texts:
I've hit an issue with using the Dialog Preserve command. I get the error message "Invalid Dialog preserve command. There is no active dialog". As far as I can tell I've followed the guidance provided in the MapBasic reference guide. Here is a subset of my code which reproduces the error in MapInfo 11.5 and 12

' * Include mapbasic.def - required for various MapBasic commands including "CMD_INFO_XXX"

Include "mapbasic.def"

' * OK button handler for Dialog

Declare Sub gmeu_goToGridReference ( gmeu_userEnteredGridRef As String )

Sub gmeu_goToGridReference ( gmeu_userEnteredGridRef As String )

    Dialog Preserve
    Print "You entered: " + gmeu_userEnteredGridRef

End Sub

Declare Sub Main

Sub Main

    ' * Create a custom dialog box to obtain the grid reference from the user

    Dim gmeu_userEnteredGridRef As String

    Dialog
        Title "Go to Grid Reference"
        Control StaticText
            Title "Enter a grid reference:"
        Control EditText
            Into gmeu_userEnteredGridRef
        Control OKButton
       Control CancelButton
    If CommandInfo( CMD_INFO_DLG_OK ) Then
        Call gmeu_goToGridReference( gmeu_userEnteredGridRef )
    End If

End Sub

I've tried moving the position of the declare statements from the top to just before the Subs are coded but this has not changed the result - I didn't think it would but wanted to be thorough.


texts: Errors importing shapely.geometry into QGIS 3.6

texts:
I got error messages when executing 2nd statement below in python console of QGIS 3.6 (and QGIS 3.8)

import shapely

from shapely.geometry import Polygon

If I execute the 2nd statement again, there would be no error messages. After the 2nd execution, all the functionalities (Polygon, LineString, ...) used in my python scripts appear working properly.

Is there any way I can avoid the error messages?
(My QGIS 2.18 arrears no problems with shapely.)

The first and last few lines of the error messages are below:

--- Logging error ---
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.6\apps\Python37\lib\logging\__init__.py", line 985, in emit
    stream.write(msg)
AttributeError: 'NoneType' object has no attribute 'write'
Call stack:
  File "C:/PROGRA~1/QGIS3~1.6/apps/qgis/./python\console\console.py", line 579, in runScriptEditor

The last few lines are:

File "C:\PROGRA~1\QGIS3~1.6\apps\Python37\lib\site-packages\shapely\geos.py", line 130, in <module>
    os.path.join(sys.prefix, "Library", "lib", "geos_c.dll"),
  File "C:\PROGRA~1\QGIS3~1.6\apps\Python37\lib\site-packages\shapely\geos.py", line 50, in load_dll
    LOG.debug("DLL: %r", dll)
Message: 'DLL: %r'
Arguments: (<CDLL 'C:\PROGRA~1\QGIS3~1.6\bin\geos_c.dll', handle 7fed56b0000 at 0x17fa9160>,)


texts: Activating GRASS utilities in QGIS Processing toolbox via text descriptors

texts:
This is a new follow-on question from this one. Following the documentation outlined in the answer, it is fairly easy to set up a descriptor text file for standard utilities that exist but aren't in the Processing tool box.  I have been testing this with a few official but undescribed utilities.

On restarting QGIS, the new descriptor is recognised and a dialog box will appear as expected. The process will run and, if the Options and Flags are specified in the descriptor as per the utility's manual page, the process will run without any reported error but the output is empty (in one example - r.out.bin).  In another example, despite specifying the output file I get an error stating that no output file was specified (r.out.gdal)!  This latter error I have noticed in another official utility with an official descriptor (r.out.vrml).

So, are there just some GRASS utilities that cannot run via QGIS, or am I doing something fundamentally wrong (if it's the latter then I don't understand why the official descriptors sometimes seem 'broken')?  I am using QGIS 2.18 and have tested with both GRASS 6.4.3 and 7.2 on Windows 10 (64 bit).

Here is my descriptor for r.out.bin:

r.out.bin
Exports GRASS raster map to a binary array
Raster (r.*)
ParameterRaster|input|input|False
ParameterBoolean|-i|Generate integer output|True
ParameterBoolean|-f|Generate floating-point output|False
ParameterBoolean|-h|Export array with GMT compatible header|False
ParameterBoolean|-b|Generate BIL world and header files|True
ParameterBoolean|-s|Byte swap output|False
ParameterNumber|null|Value to write out for null|None|None|0
ParameterSelection|bytes|Number of bytes per cell|1;2;4;8|0
ParameterSelection|order|Output byte order|big;little;native;swap|2
OutputFile|output|Output

EDIT:
After more research, I am beginning to think that r.out.bin itself is broken.  I have tested with several different datasets in GRASS 7.2 (standalone - not via QGIS) and r.out.bin consistently returns the same error "Unable to create file ".  This is odd as r.out.bin is apparently based on r.out.ascii and that utility works just fine. Can anybody confirm that r.out.bin actually works?

EDIT2:
After further testing, r.out.bin does work OK in standalone GRASS if a non-default file name/location is chosen.  However, this does not solve the problem when running through QGIS Processing.


texts: Cannot solely remove M value of Shapefiles

texts:
I have a couple of Shapefiles which are of geometry type ZM (like PolygonZM, MultiLineStringZM, etc.), however my database (a PostGIS database) is of geometry type Z (like PolygonZ, MultiPointZ, etc.). When I import my shapefile into database, it will of course report an error. However, when I convert these Shapefiles, I find that M values cannot be removed without removing Z value.
Here are my efforts:
1. In ArcMap, use Feature Class to Feature Class in ArcToolbox, click "Environment" and disable M value there. However, output Shapefiles still have M values.
2. Using Export in QGIS, and choose PolygonZ or other geometry types with Z only in SHPT field. Again, output Shapefiles still have M values.
3. Use ogr2ogr. Here is the command line:

ogr2ogr -nlt polygonz new.shp old.shp

Again, output Shapefiles still have M values.
However, if I discard Z values as well:

ogr2ogr -nlt polygon new.shp old.shp

Output Shapefiles will have neither Z nor M values.
But when I try to add Z again, like:

ogr2ogr -nlt polygonz new1.shp new.shp

Both Z and M will be added, thus returning to PolygonZM.
Since these shapefiles have no actual Z values at all, I can safely discard Z values along with M values if I only need to deal with these files themselves. However, in my database, geometry columns are storing Z values, and they cannot be discarded. Can anyone tell me why can't I just discard M values and preserve Z values for these Shapefiles?


texts: Communicating between Tim McGee' Search and Attribute Table?

texts:
I am trying to incorporate Tim McGee's attribute table and search widgets into my existing CMV site.  This is what I have done so far:

The search tool is in my TOC
The bottom pane is configured for table
I added a query to display upon opening, just to test the table, and the results of the query DO display in the table
I added a query based on my operational layer in the searchWidget.js and I see the query in my Search tool

When I click “Search” or set up a query by shape, though, nothing is happening.  I have only edited the viewer and searchWidget js files.

Can anyone steer me in the right direction?


texts: ArcGIS: Reclassify Values as Parameter in ModelBuilder

texts:
Is there a way in ArGIS' Model Builder to set the values for reclassification (Tool Reclassify) as Model Parameters? The Reclassification is an important part of my model and I need to be able to change the threshold values for reclassification in an easy way in the model window.

The only workaround I can so far think of, would be to create several rasters by using the Raster Calculator, using Double variables which can be set as Parameters...


texts: Converting Multiple Polylines to Polygons

texts: I have a shapefile for Parcels but the whole thing seems to be divided into polylines rather than polygons. I'm trying to create an existing land use map. I was hoping there is an easier way to convert multiple polylines to polygons rather than going through each one by one and merging them in editor?


texts: How to find a old business that doesnt exist anymore on Google Maps / Google Earth

texts: Is there a way to "travel in time" on Google Maps or Google Earth to see the businesses from a particular year from a particular location?


texts: Defining circular area around fixed location?

texts:
I was trying to define a circular area around a fixed location, with known latitude and longitude. The radius of the circle will be less than 90 m. I want to get an array of all the locations present in that area. My main aim is when two or more circles intersect then to get all the location coordinates(i.e. latitude and longitude), which are present in both the circles(i.e. the intersecting area).

That's my approach, I am a beginner in this GIS field.


texts: Route Analyst acts odd at using custom impedance

texts:
I need to do a route analysis with a custom made impedance, which is impedance based on the surface type. Here is the problem: while the route Analysis acts perfectly fine with the default cost "Length" (image 1), it will be messed up with my unknown custom cost "Surface_tpye" (image 2), which is based on the field suface_cost. Since I am using the same stops on both scenarios, there can't be breaks.

Any ideas to solve this problem?

I also tried with the field length as impedance, same problems. Only works with default cost.


texts: Where to find a Placename ID for Tennengebirge for use in GML?

texts: I am new to GML and trying to create an Instance file - I need the ID for Tennengebirge but I am not sure how to find the Inspire ID for this Protected Site?


texts: Searching 'select by color' for raster or map data

texts:
I'm planing to create an online map component with which the user can create/draw polygons based on underlying raster data.

Is there some library / tool / procedure existing which enables the creation of polygons based on the colors of the raster data (or map data) - something like the 'select by color' in photoshop?
Example: the user clicks in the center of a lake and the tool creates a polygon which covers the complete lake at the coastline.


texts:
I'm trying to install a module in GRASS GIS, the v.lidar.mcc, however the script g.extension is showing the following error below:

*

  g.extension extension=v.lidar.mcc url=
  Traceback (most recent call last):   File
  "/usr/lib/grass73/scripts/g.extension", line 1729, in
      sys.exit(main())   File "/usr/lib/grass73/scripts/g.extension", line 1704, in main
      check_dirs()   File "/usr/lib/grass73/scripts/g.extension", line 1381, in check_dirs
      check_style_files('grass_logo.png')   File "/usr/lib/grass73/scripts/g.extension", line 1355, in
  check_style_files
      shutil.copyfile(dist_file, addons_file)   File "/usr/lib/python2.7/shutil.py", line 82, in copyfile
      with open(src, 'rb') as fsrc: IOError: [Errno 2] No such file or directory: '/usr/lib/grass73/docs/html/grass_logo.png'

*

I had installed this module normally in GRASS on Windows, however the UBUNTU am unsuccessful.

Anyone know how to solve it? I tried to update GRASS, install and uninstall, install via python but the error keeps popping up!


texts: Tobler's hiking function: resulting walking speed in KmH or m/s?

texts:
I was reading some literature dealing with least-cost path calculation and accumulate cost-surface generation using the Tobler's hiking function. I was reading this literature since I am trying to implement in R the generation of isochrones around a starting location (Integrating use of surface distance when calculating cumulative cost-surface based on walking pace?).

I can't understand what seems to me quite odd: some sources (e.g., Wikipedia) indicate that the speed predicted by the function is in KmH, while the author of the gdistance R package (see this PDF at page 14) indicates that the off-path speed would be in meter-per-second (he actually indicates m/s).

I am having hard time in understanding which is correct. This is not a trivial issue since different results are bound to crop out using Kmh vs. m/s.


texts:
I want to create pansharpened image from landsat 8 composite image (R,G,B and NIR bands) in ArcGIS using Esri method in "Create Pan-sharpened Raster Dataset" tool.
For weighting I used recommended value from ESRI resource website:

Band 1 weight: 0.166
Band 2 weight: 0.160
Band 3 weight: 0.166
Band 4 weight: 0.3

I have got the result where red, green and blue bands are sharpened but NOT NIR one.

Any suggestion what I might doing wrong?


texts:
I'm using QGIS 3.16. I noticed that Apple Maps has satellite images that I would like to use as a basemap.
Can Apple Maps be added to QGIS?
Apple Maps on DuckDuckGo
Apple Maps on Satellites Pro


texts: How to identify overlapping GeoJSON layers in Leaflet?

texts:
I have a Leaflet map containing overlapping GeoJSON layers. How can I identify all features at the location clicked?

See http://jsbin.com/zavuda/edit?html,js,output for a demo using Leaflet sample data.

Note that the Campus layer can be clicked when it's the top layer, and its attributes are shown via the layer's onEachFeature function.

Then use the toggle under the map to switch on the States layer, which overlaps the Campus layer. After this the Campus layer can no longer be clicked on.

Is it possible to run an "identify" query at the location clicked by the user, and return all intersecting features?


texts: As a GIS Coordinator I am constantly amazed at how difficult it is to comprehend the literature produced by the European Commission regarding the Inspire directive. If this is hard work for people who work full time in the profession how can anyone expect anyone working on the peripheries of GIS to understand it´s demands? With this in mind, does anyone know of any literature describing INSPIRE for the common man? I saw an "INSPIRE for dummies" pdf once which had great reviews but it was unfortunately in french. If anyone has something similar from a good source it would be very helpful. I am of the opinion that INSPIRE doesn´t have to be difficult. It has just become too complicated through bureaucracy.


texts: QGIS 2.18.15 crashes when printing to pdf from Mac?

texts:
I'm a beginner in QGIS.

I've tried printing a vector project (not too large) to pdf but the program crashes every time. Here is some of the information provided by the crash report generated by my Mac:

Exception Type:        EXC_BAD_ACCESS (SIGSEGV)
Exception Codes:       KERN_INVALID_ADDRESS at 0x0000600081c2036d
Exception Note:        EXC_CORPSE_NOTIFY

Termination Signal:    Segmentation fault: 11
Termination Reason:    Namespace SIGNAL, Code 0xb
Terminating Process:   exc handler [0]

I wasn't able to find an answer that seems to apply to my issue in the previous questions.


texts:
I have 4 shapefiles which are 2 locations, roads, and bridges.
I am trying to find the best way from A to B with fewer bridges.

I tried network analysis but I don't have a Network data source, should I create one or there is a different way to solve it?


texts: What's the benefit of heterogeneous geometry column in PostGIS?

texts:
I read that in PostGIS, one can have different types of geometry, say, POINT and LINESTRING?? in the same column/field of a table.

I have used ESRI products a lot where each table can only contain one type of geometry. I am just curious what is benefit/motivation for allowing mixed geometry types in PostGIS columns?

To me, what PostGIS did is to have an array of Java ObjectS, which is kind of pointless, since an Object (or Geometry in the case of PostGIS) can be anything. We can't really reason much about the base class geometry (e.g. its dimension and inoperability with other entities etc.) because it's too abstract.

For example, polygons have no length, points have no meaningful lengths or areas (because points are often simplified representation of linear/areal objects). Naturally, one can't do statistics/aggregates such as finding out the average measure of a geometry field because a measure may not be defined for some rows, due to the lack of constraint on the specific type of the geometry.

Am I missing some beneficial aspects of the heterogeneous columns?

Also, are there any practical examples of hetergeneous columns for non-geometric types in Postgres? For example, a column that can be both float-ing point numbers and bytea, or is the mixed columns a PostGIS invention?


texts: how do you do a buffer analysis of a geographical location?

texts: I want to use a buffer analysis to identify the geodemographics of Incinerators. I have downloaded the Output Area Classification (OAC) map and also know the locations of all the Incinerators I want to analyse from an online map. I want to know how i can perform a buffer analysis from the geographical point in which the Incinerators are located based on the OAC map.


texts: How to compare accuracy of two DEM using GRASS

texts: How do I compare accuracy of two DEM data (lidar and srtm) (possibly  through rmse) using GRASS GIS?


texts:
I want to know how I can show that write on my search box and appear in a rectangle inside map (sync up two elements, search box and the Text area).

Like this:

[link]
(https://design.ateljekartan.com/?lat=40.4167754&lng=-3.7037902&zoom=15&text=Madrid&orientation=landscape&styleName=fullscreen-adriatic)


texts: Do I need to transform wind vectors when changing projection from EPSG4326 to EPSG24373?

texts:
I have a set of points and polygons given in "normal" lat/lon coordinates EPSG4326, which I need to transform to EPSG24373.  I can do this without problem using PROJ4.  However, I also have wind vectors (given as u and v, i.e., m s-1 in N and E direction) for several point locations.  Question is:  When I go from EPSG4326 to EPSG24373, do I also need to transform the wind vectors in any way?

My guess would be no, because the wind is given in m s-1 anyways (which is basically the target projection, except for the difference Indian yard <-> meter).  However, I wanted to have an expert opinion on this, as I'm not too familiar with geographical projections.


texts:
I have a number of bus stops (green points) and a layer of population data (square polygons) containing information on the number of households.

I want to know how many households are within the 500m buffer of the points.

How do I sum the attribute (number of households) to a new value for each point?


texts:
I'm trying to display some GeoJSON files with D3. The code is pretty simple:

import europe_data from "./data/europe.json";
import areas_data from "./data/LV.json";

/******************************************/

const [states, setStates] = useState(europe_data);
const [areas, setAreas] = useState(areas_data);

/******************************************/

const projection = geoMercator().fitSize([width, height], states);
const pathGenerator = geoPath().projection(projection);

// render each state
svg
  .selectAll(".state")
  .data(states.features)
  .join("path")
  .attr("class", "state")
  .attr("d", (feature) => pathGenerator(feature));

// render LV areas
svg
  .selectAll(".area")
  .data(areas.features)
  .join("path")
  .attr("class", "area")
  .attr("d", (feature) => pathGenerator(feature));

However, I can't correctly display LV.

This is the rendering of Europe.

This should be the rendering of some areas in Latvia, filled in purple.

I need to work with files from various sources and format, so I'm converting all of them with mapshaper in GeoJSON file formatted in WGS84.

These are the GeoJSON files I'm using.

Why can't I display LV.json?

EDIT:

The GeoJSON files both contain the keys type, geometry, and properties, even if in a different order. It doesn't affect the output.


texts: Extract Administrative Boundaries into OSM poly Format

texts:
I need to extract some administrative boundaries from big OSM extract, and put them in .poly file.
Are there any tools for this task?

I've tried Perl script by Gary68 (http://svn.openstreetmap.org/applications/utils/gary68/boundaries.pl), but it returned something like
ERROR: at least one needed node missing in osm file.
ERROR: at least one needed way missing in osm file.
ERROR: at least one needed node or way missing.

Also I've installed osmium, but I don't see tool for my task.

What would you advise me?

(Basically I need to get a list of all houses on all streets of a certain city, and do it automatically for number of cities)


texts: Add multiple fields to a batch of layers in QGIS?

texts:
I have multiple polygon .shp files that are denoting farm boundaries. The sources are varied and the original fields vary in name and number. These boundary files are to be imported into a Farm Management Software and require specific fields for identifying the data on import. The data in the fields will have to be entered manually. I am using QGIS 2.18

I have limited knowledge of scripting. A plugin would be awesome. I am aware of the "Add Field" in the Processing Toolbox. I could add each field separately through the batch function.

Is it possible to add more than one field to a layer in a batch configuration?


texts: Clip Spatial object to bounding box in R

texts:
Given a Spatial object in R, how do I clip all of its elements to lie within a bounding box?

There are two things I'd like to do (ideally I'd know how to do both, but either is an acceptable solution to my current problem--restricting a polygon shapefile to the continental U.S.).

Drop each element not fully within the bounding box.  This seems like bbox()<- would be the logical way, but no such method exists.
Do a true clip operation, such that non-infinitesimal elements (e.g. polygons, lines) are cut off at the boundary.  sp::bbox lacks an assignment method, so the only way I've come up with would be to use over or gContains/gCrosses in conjunction with a SpatialPolygons object containing a box with the new bounding box's coordinates.  Then when clipping a polygon object, you'd have to figure out which are contained vs. cross, and alter the coordinates of those polygons so that they don't exceed the box.  Or something like gIntersection.  But surely there's a simpler way?

While I know that there are many problems with bounding boxes, and that a spatial overlay to a polygon that defines the region of interest is typically preferable, in many situations, bounding boxes work fine and are simpler.


texts: Performing spatial join between lines that travel along each other in ArcGIS Desktop?

texts: I have 2 layers, one representing streets, and one representing pipes. The pipe layer has a street field that I want it to be automatically filled (I want the pipe to show the street that it's on)  I tried performing a spatial join based on 'closest' but I encountered the following problem: many pipes never intersect the center of the streets they're on (they are usually parallel), and, at the spatial join, those respective pipes don't get the needed data from the street they go along but from the ones that are perpendicular (cause they intersect the pipes, so they have one point that is closest).


texts: How do a transform a geotiff image to use the projection of another image?

texts:
I have two sets of glib2 data from different sources which I would like to process together.

Each image has a different projection and resolution (I don't know if resolution is the correct word to use here.)

for example, here are my two images:

and

both of these images I generated by parsing the glib2 data and writing the output to geotiff. Here's my code:

data = gdal.Open(file)
driver = gdal.GetDriverByName('GTiff')

image = driver.Create(output_file, sizex, sizey, 4, gdal.GDT_Byte)

image.SetGeoTransform(data.GetGeoTransform())
image.SetProjection(data.GetProjection())

band = data.GetRasterBand(1).ReadAsArray()
shape = numpy.shape(band)

r = image.GetRasterBand(1)
g = image.GetRasterBand(2)
b = image.GetRasterBand(3)
a = image.GetRasterBand(4)

aarray = a.ReadAsArray()

#.. scale arrays

r.WriteArray(band)
r.FlushCache()
g.WriteArray(band)
g.FlushCache()
b.WriteArray(band)
b.FlushCache()
a.WriteArray(aarray)
a.FlushCache()

image = None

I have tried parsing in both sets of glib2 data and then at the point where I'm setting the projection and geotransform I use the other set of transforms than the data.

Like this:

data_1 = gdal.Open(file_1)
data_2 = gdal.Open(file_2)

image.SetGeoTransform(data_2.GetGeoTransform())
image.SetProjection(data_2.GetProjection())

band = data_1.GetRasterBand(1).ReadAsArray()

So the data used is different from the transform and projection used.

When I do this, it makes no visible difference. It seems like the change in the transform and projection have no impact.

I know I'm way off base here, because I think I need to use the same transform and projection data as I do the actual data, but I don't know how to transform the image to another projection.

How do I do this?


texts: Import error when utilizing QGIS's Python Init function

texts:
Based on the answer to Referencing Qt UIs external to plugin with python. I'm trying to utilize the "Layer Properties / Fields / Python Init function" to execute a file (test_file.py) that is located in the folder of the plugin that I'm making.

I've made sure that the path to my plugin folder is in PYTHONPATH via "Settings / Options... / System / Environment" ... and restarted QGIS.

When the layer, for which I'm utilizing Init function, is initialized and a feature is selected test.py is executed. However, I get three consecutive error messages saying:

1) An error occured during execution of following code:
      import test_file
      Traceback (most recent call last):
        File "", line 1, in
        File "C:/PROGRA~1/QGISVA~1/apps/qgis/./python\qgis\utils.py", line 453, in _import
          mod = _builtin_import(name, globals, locals, fromlist, level)
      ImportError: No module named test_file

2) An error occured during execution of following code:
      if hasattr(test_file,'DEBUGMODE') and test_file.DEBUGMODE:    reload(test_file)
      Traceback (most recent call last):
        File "", line 1, in
      NameError: name 'test_file' is not defined

3) An error occured during execution of following code:
      _qgis_feature_form_20140526104937007 = test_file.py(_qgis_featureform_2, _qgis_layer_workgroups20140526104723778, _qgis_feature_20140526104937007)
      Traceback (most recent call last):
        File "", line 1, in
      NameError: name 'test_file' is not defined

Any suggestions how to fix this?


texts: OpenLayers - WFS Projection not showing EPSG:25832

texts:
I can't solve the display of my WFS.
I can see the request returning the data in my debugger, but it is not showing.
My reference for solving this problem was Openlayers3 using custom EPSG:25832.

Am I right, suspecting that I do not need to further transform my vectorLayer?
Is my center causing the problem, because I transform to 'EPSG:3857' rather then 'EPSG:25832' (again I would suspect not - since it is just the center)?

let myProjectionName = 'EPSG:25832';
proj4.defs(myProjectionName, '+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs');
let myProjection = ol.proj.get(myProjectionName);
this.mapView = new ol.View({
  center: ol.proj.transform([9.98, 53.55], 'EPSG:4326', 'EPSG:3857'),
  zoom: 9,
  projection: myProjection
});

let vectorLayer = new ol.source.Vector({
  format: new ol.format.GML(),
  url: function(extent) {
    return 'https://geodienste.hamburg.de/HH_WFS_Statistik_Stadtteile_Wahlergebnisse' +
      '?version=1.1.0&request=GetFeature&typename=Statistik_Stadtteile_Wahlergebnisse:Buergerschaftswahl_15.02.2015_-_Wahlbeteiligung_in_Prozent' +
      '&bbox' + extent.join(',') + ',EPSG:25832';
  },
  strategy: ol.loadingstrategy.bbox,
  projection: 'EPSG:25832'
});

this.map = new ol.Map({
  target: this.mapId2,
  layers: [
    new ol.layer.Tile({
      source: new ol.source.OSM()
    }),
    new ol.layer.Vector({
      source: vectorLayer
    })
  ],
  view: this.mapView
});


texts: Highlight L.divIcon on mouseover or programmatically in Leaflet map

texts:
I want to highlight L.divIcon svg markers on mouseover and/or from an "outside action" like pressing a button.

here is a simplified testcase https://jsfiddle.net/sxvLykkt/5/

Markers are generated dynamically (geoJson originally) and added to a L.FeatureGroup().
On mouseover I set a bigger version of the icon (divIconActive) on a temporary layer, that I remove on mouseout. Unfortunately this doesn't work as indented (it's flickering on mouseover, the event is firing mouseover and mouseout I believe). How can I solve this?

And how can I access the markers when clicking one of the buttons? Somehow over their index I believe?! I can't wrap my head around.

Here is a portion of the code how the markers are created:

// init map and tileLayer -> jsfiddle
var coords = [[53, 13],[49, 10],[46, 12],[51, 16]];

$.each(coords, function(i,e){
  // create the button
  $('#controls').append('<button>'+i+'</button>')

  var marker = L.marker(e, {
    icon: divIcon,
    id: i
  });

  locationLayer.addLayer(marker);

  marker.on('mouseover', function(e){
    markerTemp = L.marker(e.latlng, {
        icon: divIconActive
    }).addTo(map);

  });

  marker.on('mouseout', function(e){
    markerTemp.remove();
  });

});

locationLayer.addTo(map);

$('button').on('click', function(e){
    alert('Highlight the right marker!')
});

PS: this is a cross post from Stackoverflow.

Edit: JsFiddle https://jsfiddle.net/sxvLykkt/12/ with the proposed, working solution by TomazicM


texts: Count number of pixels from .msk in Python to derive study area size

texts:
Apologies I am totally new to GIS in Python.

I am trying to find out the total area of the Barents Sea mask the NSIDC use to calculate their products and therefore I need to know how many pixels the mask is (so I can derive total area since lat and lon resolution = 25km). However the mask is a numpy .msk format which I am not sure how to process.

The .msk file they use can be downloaded from their FTP
"Arctic_region_mask_Meier_AnnGlaciol2007.msk"> ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/seaice_analysis/

and is loaded in

import numpy as np
regional_mask = np.fromfile('./Arctic_region_mask_Meier_AnnGlaciol2007.msk',dtype=np.uint8).reshape((448, 304))

The Barents Sea has the mask value of 8 in the above .msk file, but I do not know how to just find the total pixel count from this value?

Additionally since the .msk does not contain lat/lon I do not know how to reconstruct their lat and lon from a .dat file?

"psn25lats_v3.dat" and "psn25lons_v3.dat" > ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/tools/


texts: QGIS Extend Toolbar (GUI) Button Functionality via Python Plugin?

texts:
Is it possible to piggyback onto the QGIS split method via a python plugin?  Are toolbar buttons accessible via the python API?

This looks like the class that I need to access.  Access ideas?

Proposed Workflow to Give an Implementation Example:

Initialize editing
Split the line as per How to Split Features in QGIS?
Launch a custom GUI or otherwise call a set of functions to start to build a list that I push to the DB table on saving edits.


texts: Calculating temperature of given area using MODIS LST data?

texts:
I am using Google Earth Engine and the MODIS Land Surface Temperature data to try and calculate the temperature of a square of land.

Here is my code so far:

var dataset = ee.ImageCollection('MODIS/006/MOD11A1')
              .filter(ee.Filter.date('2018-12-01', '2018-12-10'));
var landSurfaceTemperature = dataset.select('LST_Day_1km');
var landSurfaceTemperatureVis = {
  min: 13000.0,
  max: 16500.0,
  palette: [
    '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
    '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
    '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
    'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
    'ff0000', 'de0101', 'c21301', 'a71001', '911003'
  ],
};
Map.setCenter(6.746, 46.529, 2);
Map.addLayer(
    landSurfaceTemperature, landSurfaceTemperatureVis,
    'Land Surface Temperature');

var temperature = landSurfaceTemperature * 0.02 - 273.15
print(temperature)

As you can see, I am trying to calculate the temperature using the formula T = DN * 0.02 - 273.15 (I found this formula in the following thread How do I convert the LST values on the MODIS LST Image to degree celsius
).

When run, the program produces an output of 'NaN'. I believe I am trying to use the data in the wrong format, but I do not know the correct one.


texts: Rasterio doesn't display the differnet bands of image

texts:
I have one raster which has three different bands.
The two first bands look like this (in qgis):

and the third band looks like this ( missing values in some places:

The three bands sit on top of each other.
The problem is that whenever I try todisplay those layers with  rasterio it always shows the "extent" of the third layer:

img=rasterio.open("21072019.tif")
show(img,1)

#display other bands:
show(img,2)

show(img,3)

with all those bands the image looks exactly the same, th eonly thing that changes is the axis.
My question: How can I display the different bands? why is this hapenning? do the difffeernt bands has different shape?


texts: How to calculate surface out of GPS points in QGIS

texts:
I have uploaded several squares (which consist out of 4 GPS points, not connected). I  need the total surface of the squares.

Is there a simple way of doing this?

Do I create polygons and than calculate the surface?


texts:
I'm working on a project in QGIS right now that uses Canada and Canadian waterbodies as a basemap. However, the waterbodies file I was provided seems to warp in the Northwest Territories when you zoom out to see the complete map.

I've looked for a replacement shapefile online but I haven't found a suitable replacement so I'm trying to fix the issue instead.

Apologies if I've left anything out, this is my first post.


texts: How to extract values associated with identified spatial points that is extracted from a spatial polygons data frame using buffer ins

texts:
I plotted domestic livestock abundance (each spatial point corresponds to different values, it appears as 1,2.3......but each of these 1, 2,3 has corresponding number of livestock counts in it) in a 4 by 4km grid (raster layer).

Then i plot the facal sample collection locations and around each sample made a buffer of 10km radius of spatialpolygonsdataframe. then got the list of or number of domestic livestock abundance in each buffer. But now stuck with how to extract the values associated with each of extracted lists.

extent(r)
res(r) <- 4000.0 #sets the resolution

r <- setValues(r, 1:ncell(r))
grid <- as(r,'SpatialPolygonsDataFrame') # grid for making SPDF

wp.per.grid <- over(wp.count.sp, grid) # gives count per grid

buff <-  gBuffer(scats.sp, width=10000, byid = TRUE) # made buffer

d1<-data.frame(pointID = 1:nrow(wildp.count.sp))
for(i in 1:length(buff@polygons)){
  column<-as.logical(point.in.polygon(point.x = wp.count.sp@coords[,1],
                                      point.y = wp.count.sp@coords[,2],
                                      pol.x = buff@polygons[[i]]@Polygons[[1]]@coords[,1],
                                      pol.y = buff@polygons[[i]]@Polygons[[1]]@coords[,2]))
  d1<-cbind(df,column)
  colnames(df)[i+1]<-paste0("buffer",i)
}
# gives no of counts in each buffer created
d2<-data.frame(Buffer = 1:length(buff@polygons), numberOfpoints = NA)
for(i in 1:nrow(d2)){
  d2$numberOfpoints[i]<-sum(as.integer(df[,(1+i)]))
}

how to find out which are these buffers esp for those that are more than 1 and how to get the no.of counts associated with each of the identified points inside buffer???


texts:
I'm tryng to clip a raster with a shapefile using the tool in QGIS 2.18: Raster--> Extraction--> Clip.

Since I want not only to clip the raster layer, but also clip the dataset in order to have a lighter file (using the shape as a mask and checking the box "clip also the dataset" as you can see from the image [forgive me, in the image everything is in Italian]).

I tried to solve my problem with previuos discussion but nothing. In fact, the most similar was this one Rasterfile clipping error: Cannot compute bounding box of cutline, where the solution was to set the same CRS to the raster and the shape. I checked and it was the same CRS, but the error continued to appear.

If I try to clip the layer without checking the box for clipping also the dataset, it works.

Any suggestion?


texts: Attribute data different for duplicated layer?

texts:
I have a layer that I want to style in different ways for different print composers. When I duplicate the layer it uses the parent shapefile, but changes made to the attribute tables aren't reflected across duplicate layers.

Is this expected behaviour? Do duplicated layers have separate attribute tables?


texts: Move/Drag multiple graphics in ArcGIS JavaScript API

texts:
In ArcGIS JavaScript API, I want to click a GraphicsLayer object (Graphic) on the map and perform a dragging MOVE operation for it as well as some of its neighboring graphics.  I am looking for the visual effect of several graphic objects moving together.  The EDIT toolbar MOVE operation accepts a single graphic, not multiple, as far as I see.

editToolbar.activate(Edit.MOVE, evt.graphic);

Can I aggregate into a single Graphic for use with EDIT toolbar?

Activate a toolbar for each Graphic?