texts: Fixing Raster Calculator Error 000539 ImportError: No module named numpy?

texts:
I've created 4 hillshades for a particular area to represent the 4 seasons in a year. The equation I now need to use to determine solar insulation is 1000*hillshade/255 in raster calculator. However, I keep getting this error:

  Messages
  Executing: RasterCalculator "(1000 * "Spring_HS") / 255" "E:\Geog 170\Lab 2\DEM\spring_insul"
  Start Time: Sat Apr 11 09:44:25 2015
  (1000 * Raster(r"Spring_HS")) / 255

  ERROR 000539: Error running expression: rcexec()
  Traceback (most recent call last):
    File "", line 1, in
    File "", line 2, in rcexec
    File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy__init__.py", line 24, in
      from arcpy.toolbox import *
    File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\toolbox.py", line 356, in
      from management import Graph, GraphTemplate
    File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 22, in
      import _management
    File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy_management.py", line 14, in
      import _graph
    File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy_graph.py", line 27, in
      import numpy
  ImportError: No module named numpy

Failed to execute (RasterCalculator).
Failed at Sat Apr 11 09:44:25 2015 (Elapsed Time: 0.12 seconds)

I'm not entirely sure what the error message means, but a few people, it seems, determined that it is due to an issue with the raster calculator and arcpy. Please help me to solve this issue or find some kind of work around. Additionally, as a student, it would be great if we could also try to dissect this issue so that I may better understand just what the heck is going on.

It seems that the error code being put out is referring to an issue importing something called "numPy". Here is a duplicate of this question with a slew of different solutions than the solution found in this thread.


texts: r - sf::st_within output is sparse geometry binary predicament format. How to convert to character or logical vector

texts:
The same question is asked here:
https://stackoverflow.com/questions/49294933/r-convert-output-from-sfst-within-to-vector
but I thought this might be a more suitable place for GIS related questions.

Im trying to use the sf package in R to see if sf object is within another sf object with the st_within function. My issue is with the output of this function which is sparse geometry binary predicate - sgbp and I need a vector as an output so that I can use the dplyr package afterwards for filtering. Here is a simplified example:

# object 1: I will test if it is inside object 2
df <- data.frame(lon = c(2.5, 3, 3.5), lat = c(2.5, 3, 3.5), var = 1) %>%
st_as_sf(coords = c("lon", "lat"), dim = "XY") %>% st_set_crs(4326) %>%
  summarise(var = sum(var), do_union = F) %>% st_cast("LINESTRING")

# object 2: I will test if it contains object 1
box <- data.frame(lon = c(2, 4, 4, 2, 2), lat = c(2, 2, 4, 4,2), var = 1) %>%
  st_as_sf(coords = c("lon", "lat"), dim = "XY") %>% st_set_crs(4326) %>%
  summarise(var = sum(var), do_union = F) %>% st_cast("POLYGON")

# test 1
df$indicator <- st_within(df$geometry, box$geometry) # gives geometric binary predicate on pairs of sf sets which cannot be used
df <- df %>% filter(indicator == 1)

This gives  Error: Column indicator must be a 1d atomic vector or a list.

I tried solving this problem below:

# test 2
df$indicator <- st_within(df$geometry, box$geometry, sparse = F) %>%
  diag() # gives matrix that I convert with diag() into vector
df <- df %>% filter(indicator == FALSE)

This works, it removes the row that contains TRUE values but the process of making a matrix is very slow for my calculations since my real data contains many observations. Is there a way to make the output of st_within a character vector, or maybe a way to convert sgbp to a character vector compatible with dplyr without making a matrix?


texts:
I'm trying to associate a previously created style with a layer, according to https://docs.geoserver.org/latest/en/api/#1.0.0/styles.yaml, through a NodeJS code.
POSTing to:

http://...:8080/geoserver/rest/layers/camada_de_estacoes_3/styles?default=true

Body:
{ name: "camada_de_estacoes_3", filename: "camada_de_estacoes_3.sld" }

PS1: yes, style with same name as layer;
PS2: also tried with "[workspace name]:" in front of names of things, no success;
PS3: Geoserver version 2.13.2.
Response: status 500 - "name"
What am I doing wrong? Thanks in advance!
Log piece:

2020-09-29 12:00:22,490 ERROR [geoserver.rest] - name
com.thoughtworks.xstream.mapper.CannotResolveClassException: name at
com.thoughtworks.xstream.mapper.DefaultMapper.realClass(DefaultMapper.java:81)
at
com.thoughtworks.xstream.mapper.MapperWrapper.realClass(MapperWrapper.java:125)
at
com.thoughtworks.xstream.mapper.DynamicProxyMapper.realClass(DynamicProxyMapper.java:55)
at
com.thoughtworks.xstream.mapper.MapperWrapper.realClass(MapperWrapper.java:125)
at
com.thoughtworks.xstream.mapper.PackageAliasingMapper.realClass(PackageAliasingMapper.java:88)
at
com.thoughtworks.xstream.mapper.MapperWrapper.realClass(MapperWrapper.java:125)
at
com.thoughtworks.xstream.mapper.ClassAliasingMapper.realClass(ClassAliasingMapper.java:79)
at
com.thoughtworks.xstream.mapper.MapperWrapper.realClass(MapperWrapper.java:125)
at
com.thoughtworks.xstream.mapper.ArrayMapper.realClass(ArrayMapper.java:74)
at
com.thoughtworks.xstream.mapper.MapperWrapper.realClass(MapperWrapper.java:125)
at
com.thoughtworks.xstream.mapper.SecurityMapper.realClass(SecurityMapper.java:71)
at
com.thoughtworks.xstream.mapper.MapperWrapper.realClass(MapperWrapper.java:125)
at
org.geoserver.config.util.SecureXStream$DetailedSecurityExceptionWrapper.realClass(SecureXStream.java:175)
at
com.thoughtworks.xstream.mapper.MapperWrapper.realClass(MapperWrapper.java:125)
at
com.thoughtworks.xstream.mapper.CachingMapper.realClass(CachingMapper.java:47)
at
com.thoughtworks.xstream.core.util.HierarchicalStreams.readClassType(HierarchicalStreams.java:29)
at
com.thoughtworks.xstream.core.TreeUnmarshaller.start(TreeUnmarshaller.java:133)
at
com.thoughtworks.xstream.core.AbstractTreeMarshallingStrategy.unmarshal(AbstractTreeMarshallingStrategy.java:32)
at com.thoughtworks.xstream.XStream.unmarshal(XStream.java:1486) at
com.thoughtworks.xstream.XStream.unmarshal(XStream.java:1466) at
com.thoughtworks.xstream.XStream.fromXML(XStream.java:1346) at
org.geoserver.config.util.XStreamPersister.load(XStreamPersister.java:667)
at
org.geoserver.rest.converters.XStreamJSONMessageConverter.readInternal(XStreamJSONMessageConverter.java:59)

(...)


texts: arcpy.Select_management invalid sql expression - VDI

texts:
I am using ArcMap 10.2.2 and trying to make separate feature classes from a larger feature class using a field delimiter and SQL statement but am getting an error that the SQL statement is invalid. Here is the snippet of code I'm having a problem with:

# This script looks at a featureclass(FC)/service location(SL)
# and splits it by day for web service data preparation.
# In this manner, end users can toggle off and on days of service.

#import modules
import arcpy
import os

#create input FC/SL -- the mega FC/SL
inputFC = raw_input("Paste Route FC from ArcCatalog: ")

#set default geodatabase
arcpy.env.workspace = raw_input("Paste default file GDB from ArcCatalog: ")

#designate output location as a variable so it can be used in system path
#output filename
outLocation = arcpy.env.workspace

# Designate the type of fleet, service, commodity etc
fleetType = raw_input("Name of output (e.g. Residential_Garbage, NO SPACES!!!): ")

#make a copy of the input FC so nothing gets corrupted
arcpy.MakeFeatureLayer_management(inputFC, "Copy_inputFC")

dayList = []
while True:
    day = raw_input("Type field delimiter and press enter. Type and enter 'Done' when finished: ")
    if day.lower() == 'done':
        break
    dayList.append(day)

#Make multiple FCs based on route number
for day in dayList:
    #Add field delimitters
    delimFld = arcpy.AddFieldDelimiters("Copy_inputFC", "Exist_Partition")
    print "step1"
    if day == '1':
        dayName = 'Monday'
    elif day == '2':
        dayName = 'Tuesday'
    elif day == '3':
        dayName = 'Wednesday'
    elif day == '4':
        dayName = 'Thursday'
    else:
        dayName = 'Friday'

    #Get output name
    outName = "{0}{1}{2}".format(fleetType,"_",dayName)
    outFC = os.path.join(outLocation, outName)
    print "step2", outFC

    #If routes are numbers:
    try:
        sql = "{0}={1}".format(delimFld, day)
        arcpy.Select_analysis("Copy_inputFC", outFC, sql)

    #If routes are strings
    except:
        sql = "{0} = '{1}'".format(delimFld, day)
        arcpy.Select_analysis("Copy_inputFC", outFC, sql)

    print outName

print "Finished"

This is what the error is saying:

ExecuteError: ERROR 000358: Invalid expression "Exist_Partition" = '1'
Failed to execute (Select).

However, if I do the SQL statement in ArcMap it appears to be fine:


texts: Save QGIS canvas as image with world file by python code

texts:
I have to make georeferenced TIF files from some WMS layers.
Using built-in Project SaveAsImage makes file, that is too complicated for my fieldwork program Topos. My program can't read rasters with 4 channels and 32 bit, like QGIS makes.
I want to use something like

img=QImage(QSize(800,600), QImage.Format_RGB16)

Or any other solution to change raster options.
Or any solution to make some plugin for this work.


texts: Why does setSubSymbol causes QGIS to crash?

texts:
I am having problems creating a multilayer line symbol with this code :

symbol = QgsLineSymbolV2.createSimple({})
symbol.deleteSymbolLayer(0)

symbolLine = QgsSimpleLineSymbolLayerV2()
symbol.appendSymbolLayer(symbolLine)

symbolArrow = QgsMarkerSymbolV2.createSimple({})
symbolArrow.deleteSymbolLayer(0)

symbolArrowLayer = QgsSimpleMarkerSymbolLayerV2()
symbolArrow.appendSymbolLayer(symbolArrowLayer)

symbolLineArrow = QgsMarkerLineSymbolLayerV2()
symbolLineArrow.setSubSymbol(symbolArrow) # Causes crash !!!

symbol.appendSymbolLayer(symbolLineArrow)

Is there something wrong in this code ?
Has anybody already had this issue ?


texts: Changing default location of project in qgis-web-client

texts:
we're using qgis-web-client to present our maps publicly over the net. I haven't worked with it a lot so I'm still playing around with it a bit. One thing I do not understand, yet: When qgis-web-client opens I just see a blank screen because the extents are not set correctly.

But how do I change them?

I also set var useGetProjectSettings = true; so actually qgis-web-client should zoom to the correct extents by default.

It's quite annoying to always click on the icon in the middle of the navigation arrows and then zoom to the correct location.


texts: Why do we need Mapserver/Geoserver to present data from spatial database to web browser?

texts:
I was trying "How I set up GeoServer on an old machine" and By @iant I was setting up Tomcat and all.

I suddenly had a question ... why do we need MapServer and GeoServer?

Can these servers pull data from the database layer as the user asks and only presenting what user has asked from the database?


texts: Restoring PostGIS 2.1. dump in PostGIS 2.2

texts:
I tryed to restore my old PostGIS dumps (Postgis 2.1, PostgreSQL 9.1) at my new system (PostGIS 2.2 and PostgreSQL 9.5)

I was using pgadmin3 to dump the databases and was now trying to restore it with pgadmin3 but got this error:

pg_restore: [archiver (db)] could not execute query: ERROR:  could not access file "$libdir/postgis-2.1": No such file or directory

Do I have to downgrade my PostGIS-Version back to 2.1? Or is there a more elegant way?

Im working on a Ubuntu 16.04 system.


texts: Python Toolbox (.pyt) Tool Output Layer Name

texts:
How can I manipulate the TOC display name of my output parameters when using a Python toolbox (.pyt) Tool class.  I tried adding a symbology layer, however this has no affect on the TOC name.  Is it possible to reference the layer(s) added as output from the execute method of a tool.  The expected results are added to map and TOC, but I can't figure out how to manipulate the TOC display name. I could loop through the layers in the TOC after exectution and check each layers datasource, but this seems problematic. I've searched the help and applied the extent of my Google Foo.

In my tool class I define my parameters as follows:

def getParameterInfo(self):
    extent = arcpy.Parameter(
            displayName='Analysis Extent', name='extent_features'
            , datatype='GPFeatureLayer', parameterType='Required'
            , direction='Input')
    extent.filter.list = ['Polygon', ]

    tilelayer = arcpy.Parameter(
            displayName='Analysis Tiles', name='tile_layer'
            , datatype='GPFeatureLayer', parameterType='Required'
            , direction='Output')
    tilelayer.symbology = os.path.join(os.path.dirname(__file__), 'pyltk/analysis_tiles.lyr')

Follow Up:
This normally would not be a problem as the output would reflect the name of the feature class.  However, I've been evaluating the use of the new SQLite workspace for integration with external libraries, etc. ArcGIS prepends the internal SQLite database name 'main' to all my feature classes and tables.  When ArcPy adds the results to the map it only adds the database name, and not the whole name.  For instance if the output feature class is called "..\analysis_tiles" ArcGIS will refer to this as "..\main.analysis_tiles", but the resulting TOC name is just "main".


texts: Use field calculator to change values from a,b,c etc. to 1,2,3 etc.?

texts:
I have a vector shapefile with a column with string values. How can I convert this to a column of integer values in another column but keep the categorization.

e.g. from:

col1,newcol2
A,
A,
B,
A,
C,

to:

col1,newcol2
A,1
A,1
B,2
A,1
C,3

I'm sure it is in the 'replace' or 'case when' syntax somewhere but I cannot find a clear method via Google etc.


texts: How to load a raster layer at the top of the QGIS ToC?

texts:
The following code inserts my new_layer just above the active layer in the QGIS ToC.

new_layer = iface.addRasterLayer(...)

Instead, I'd like my new_layer to be inserted at the top of the QGIS ToC.

Do you know how can I do it?


texts: No popup on map if no feature is selected

texts:
I've got a webmap showing data coming from a Geoserver instance. With better.WMS.js I made it possible to display popups.

Everything works nice except that the map shows popups too where there is no feature selected.

The code is hosted on GitHub. I'll try to put a working map online though one may see the problem clearly. The map shouldn't display a popup window when no feature is selected.

Does anyone has an idea how to fix that?


texts: Segmentation boundaries image in eCognition

texts:
I need to save the segmentation result (boundaries image) as a raster file using eCognition.

Anybody can give me some directions about this?


texts:
I would like to use interpolate point data to create a raster in QGIS. I am doing batch analysis so I would like to use a python script.

I found this question : how-to-add-an-interpolation-raster-from-the-python-console-in-qgis and related "answer" QGIS user manual.

I thus tried to get the QgsInterpolator function to work. But I could not find how to use it, nor is it present in the processing.alglist().

This question how-to-call-the-interpolation-plugin-from-the-python-console also provide some elements, but I could not find which parameters to enter for the RasterInterpolator.


texts: Math algorithm for n amount points in objective-c

texts:
I need some information about the Answer on an earlier Question.

I need this calculation for an iOS App. I have three  points with different Circles.
I understand the problem and your solution. But i have no idea how the function self works.
Do you have an idea which libary can help me to solve this problem (C C++)
I also looked at google. But i'm not really sure which lib can do this.

I need this for objective c. Also google had no solution for me.

Maybe you have me some Resources or a Mathematic description for me.
It would help me really when you have some Informations for me.


texts: Address to Zip+4 ArcGIS Address Locator or batch service? Free?

texts:
I have 200,000 addresses in MD for which I need zip+4.

I am searching for an address locator for ArcGIS or a (free?) batch service that will  output 9-digit zip (Zip+4) from address.  I have already geocoded these addresses so I have lat/long from street or point address for a good portion.

I do have Esri Street Map Premium service and that allows for 9 digit zip to Lat/long but not Address to Zip+4.  Odd. There is a "PostalExt" field in the output but i does not populate. very annoying.  Reference data must be available for the boundaries (or at least centroid) for this data but it is not available for use in the street map premium FGDB.

Since I'm only interested in one state (MD), my plan was to reverse geocode all permutations of possible 9 digit zips (about 6 million) to get their lat/longs, but it seems quite dirty to do a spatial join on these points to the previously geocoded 200k address points.

Any other suggestions very welcome. I have seen that USPS does not have or make available polygons for zip+4 but you can purchase tables of the address ranges per +4. Not sure if that would even be that helpful for me.

Thanks :)


texts: Creating polygon with exact measurements in QGIS 3

texts:
I want to create a square in QGIS 3 with dimensions of 50m by 50m exactly.

How can I implement that?


texts: How long should it take to load a point layer with readOGR?

texts:
I am afraid I am quite new to spatial analysis at anything beyond a superficial level, and I have only recently transitioned to utilizing R.  I have been quite pleased with the versatility of areal analysis, but I have just attempted to load my first point layer.
In a nutshell, it's taking much, much longer than anticipated.  I first tried to use readShapeSpatial (from maptools I believe), but I read somewhere that readOGR (from rgdal) was faster.  I let readShapeSpatial run for over an hour before killing it, and now readOGR has been running for at least 20 minutes.
Here is the line:
c<-readOGR('/home/.../DCSpatialAnalysis','CamaRes.shp')

The data attribute houses a frame containing character and numeric data.  There are 106,046 records and 46 columns.  The file size is only ~75MB (including the all elements of the shapefile).
I do not believe my hardware is underpowered for this task (though I am always happy with more horsepower):

OS: Ubuntu 13.04 (VMWare) on Windows 7
Make/Model: Dell/OptiPlex 790
Processor: i3 (3.3 GHz)
RAM:  16 GB installed

I should mention that I have been able to successfully parse and display an XML document from within my IPython Notebook session that resides in the same directory.  That suggests, to me, that reading from a mounted directory is not an issue.
Have I done something incorrectly here?  R seems to just hang in the "Busy" status, both in RStudio and IPython Notebook (via %rmagic).  How long should such a read in take?


texts: Make new Tilemill map match default Mapbox label styles

texts:
I want to make a map in tilemill that matches the way that labels are displayed in this example, down to the countries that get labels at each zoom level, which country gets abbrev. vs. full title, using lines to connect labels to smaller countries, etc.

That seems to be the default Mapbox label style (I notice that this Mozilla map acts exactly the same way), since there are no customizations there. However, when creating a new Tilemill map, there are no labels at all. What do I need to do to add labels in exactly the same way as that Mapbox example (i.e., the Mapbox default label styles)?


texts: How to make new polygon/feature from point which is NOT centroid?

texts:
I've spent all morning searching help sites and digging through the FAQs and worksheets but have failed miserably.

I am an ArcGIS novice so I hope this is not too much of a "newbie" question but I'm trying to work out how to make a new polygon/feature from a point which is NOT a centroid?

----                ----
| . | <--not this   |   |
----       this-->  .---

I have a series of points where the Easting and Northing listed represents the lower left hand corner of a 200m square. I know how to make a circular buffer which I can convert to a square with the Feature Envelope to Polygon tool but this isnt really what I want.

The data is of 3g strengths which is meant to be in grid squares. I cannot upload this file this due to NDA.

Any ideas?


texts: different results ogr2ogr command line vs osgeo import ogr with reading osm data

texts:
I am trying to read in openstreetmap data (osm.pbf) files into python. Initially, I've always used the ogr2ogr command line approach to convert the data I want from the osm.pbf to a shapefile and then read the shapefile using geopandas. However, I don't really need shapefiles so I thought it would be nice to load the data directly into python.

Unfortunately, I seem to miss ways/segments when directly reading it into python. Please see below the code I am using for both approaches:

# path to osm.pbf file
osm_path = 'C:\\planet_osm\\tanzania-latest.osm.pbf'

# import through python
driver=ogr.GetDriverByName('OSM')
data = driver.Open(osm_path)
layer = data.GetLayerByName('lines')

features=[x for x in layer]
railway = []
for feature in features:
    if feature.GetField('railway') is not None:
        osm_id = feature.GetField('osm_id')
        shapely_geo = shapely.wkt.loads(feature.geometry().ExportToWkt())
        service = feature.GetField('service')
        maxspeed = feature.GetField('maxspeed')
        highway=feature.GetField('railway')
        railway.append([osm_id,highway,service,maxspeed,shapely_geo])

rail_dir_import = gpd.GeoDataFrame(railway,columns=['osm_id','infra_type','service','maxspeed','geometry'],crs={'init': 'epsg:4326'})

# import through ogr2ogr
shape_in = 'C:\\planet_osm\\tanzania-railway.shp'

os.system("ogr2ogr -progress -overwrite -f \"ESRI Shapefile\" -sql \
               \"SELECT osm_id,railway FROM lines WHERE railway IS NOT NULL\" \
               -lco ENCODING=UTF-8 -skipfailures -nlt geometry "+shape_in+" "+osm_path)

rail_com_line = gpd.read_file(shape_in)

The python approach only gives me 563 segments, whereas the ogr2ogr command line approach (and then loading it into python) gives me 1060 segments. See below the output in visual form.

Python approach:
https://www.dropbox.com/s/kuxcaiiigdn1u58/python_import.png?dl=0

ogr2ogr command line approach:
https://www.dropbox.com/s/a854i0345e62fc4/ogr2ogr_import.png?dl=0

Does anyone know why the python approach misses data?

Update

If I use the following, I do get the right amount of railway segments:

    layer = data.ExecuteSQL("SELECT osm_id,service,maxspeed,railway FROM lines WHERE railway IS NOT NULL")

Still, that does not solve the problem above. It makes me think that there is some sort of caching problem when loading all features from a layer?


texts:
Possible Duplicate:
  Get list of coordinates for points in a layer?

Sorry for another basic question - my map is almost done but now my boss wants me to enter new outbreaks (that we don't have coordinates for) into a map that has animal population data for several provinces (one shapefile for each country).

At this point all I have for these newer outbreaks is the province that the outbreak was in - my boss has said that I can input the coordinates (roughly) for each outbreak by finding the centroid. I've used ftools for this, I've been importing data for outbreaks from excel files (using a delimited text layer) - for the outbreaks without coordinates (for which I only have data telling me what province the outbreak was in) I presume I'm going to have to make an Excel/DBF file with these coordinates. My boss says that I can get these coordinates by finding the centroids.

I've tried using "Polygon Centroids" via the ftools plugin, and have a new shapefile with points on where I presume the center of each province is.

What I'd ideally like is a list (such as in the attribute table) next to each province (which are listen in the attribute table of the shapefile I'm working with) with the coordinates for each centroid. Is there any way to get this done? The only other thing I can imagine is writing down the coordinates of each centroid of each province manually.


texts: Performance Scanning Big Rasters in GeoTIFF Format?

texts:
I am scanning large rasters in GeoTIFF format. The size is around 30GB and there is only one band. I'm using geotools to do that. I did some experiments on scanning the whole raster to find the best way to do this and realized that performance varies greatly. I don't know if these rasters were tiled. This is the code I used:

file file = new File(path);
GridCoverage2DReader  reader = new GeoTiffReaderCached(file);
GridCoverage2D coverage;
ParameterValue<OverviewPolicy> policy = AbstractGridFormat.OVERVIEW_POLICY.createValue();
policy.setValue(OverviewPolicy.IGNORE);
ParameterValue<String> gridsize = AbstractGridFormat.SUGGESTED_TILE_SIZE.createValue();
gridsize.setValue("512,512");
ParameterValue<Boolean> useJaiRead = AbstractGridFormat.USE_JAI_IMAGEREAD.createValue();
useJaiRead.setValue(true);

try {
      coverage = this.reader.read(new GeneralParameterValue[] { policy, useJaiRead});
    } catch (IOException e) {
      e.printStackTrace();
    }

RandomIter iter = RandomIterFactory.create(img, null);

//here I iterate either row by row or tile by tile. this is for reading
//values row by row
for (int i = 0; i < height; i++) {
    for (int j = 0; j < width; j++) {
        rl.iter.getSample(j, i, 0);
    }
}

In the experiments, I change between reading row by row or tile by tile and change between specifying a value for suggested tile like (512,512) or just leave it as is (the tile will be a row by default). Here is the four different ways and their timings:

With no tile suggested; reading values row by row: timing (sec) = 328.
With no tile suggested; reading values tile by tile: timing (sec) = 42334 (the worst).
With tile suggested (512,512); reading values row by row (ignoring the tiles): timing (sec) = 4295.
With tile suggested (512,512); reading values tile by tile: timing (sec) = 3409.

The best performance was for the first approach.

Why is that, and will it be true for all rasters?


texts: Get list of all map services in a Web App in ArcGIS Portal

texts: I'm in charge of moving our test Portal maps into production and looking for an easy way to get a list of all the layer sources in the web map.  In addition it would be great to be able to export a text file with the widgets, and settings for that web app.


texts: Setup for several origins for a WMTS in OpenLayers

texts:
I´m trying to setup a map with a WMTS as a layer with help form this example: https://openlayers.org/en/latest/examples/wmts.html

The layer works fine for the first three levels, but after the third, the layer get an offset.

If i look into my GetCapabilities I can see that the first three levels has the same "TopLeftCorner". But in the fourth, the value changes.

From:

<TopLeftCorner>6445536 95000</TopLeftCorner>

To:

<TopLeftCorner>6441440 95000</TopLeftCorner>

This also happens for level 5, 7 and 8, where the value changes.

As i´m using:  origin: ol.extent.getTopLeft(projectionExtent) the value is constant on 95000, 6445536, which make me belive thats why I get an offset.

Is it possible to use more than one set of x & y-values for the origin depending on the resolution? Or any other way to get around this problem? (not the GetCapabilities-example).


texts: How can I achieve high quality supervised classification?

texts:
For the supervised classification I use: QGIS, Semi-Automatic-Classification Plugin, a Landsat-8 Satallite Image and the CORINE Land Cover classes (http://uls.eionet.europa.eu/CLC2000/classes/index_html)

My goal is to perfom a high quality classification.

So far I know it is important to have many reference-areas with high correlation between actual surface and created class. Sounds easy! But it isn't! Because the entire landscape doesn't fit into just some classes.

For example: What is the "better" way to go to classify urban areas? Creating one large polygon and ignoring some green areas or creating many small/tiny polygons but therefor having less different pixels inside the polygon.

 ;

Also how do I figure out how many classes are needed? So far I decided to perform a level-1 CORINE classification, which means 5 classes.
What else is important to know in terms of a high quality classification?


texts: Adding feature layer to ArcGIS Online from REST API?

texts:
I need to add a feature layer to ArcGIS Online account from code.

i have created a developer account where if i want to create a feature layer i  have to upload a CSV file and it makes the feature layer.

  Problem is: i need to create the sane feature layer using my code.
  like is there any REST API or some API which ArcGIS provide this
  option?


texts: Highlight the selected feature with different colors

texts:
Is it possible to change the color of a selected feature in qgis 2.2.0 Valmiera? To highlight a feature I use:

aLayer.Select(feat.id())
box = vlayer.boundingBoxOfSelected()
qgis.utils.iface.mapCanvas().setExtent(box)
qgis.utils.iface.mapCanvas().refresh()

But the selected featue is colored yellow. Is it possible to define a different color for this feature?


texts: Mapbox Android SDK with Raster Tile Layer or mbTiles

texts: Is it possible to add a raster TileLayer or mbTiles file as basemap of Android SDK? I do not want to use Mapbox Maps. I need to replace it with custom map with tile format or mbTiles. After Android SDK version 0.7 Mapbox revoke the mbTileLayer option. Is there any way to add this?


texts:
I've got a Silverlight client that serializes a FeatureSet and sends to an SOE.  I'd like for the SOE to deserialize the json into an IRecordSet.   ESRI.ArcGIS.SOESupport.Conversion has a method for going from IRecordset to json, but not the reverse.

How can I deserialize json to create an IRecordset?

{
    "displayFieldName": "",
    "fieldAliases": {
        "OBJECTID": "OBJECTID",
        "DissolveID": "DissolveID",
        "CatchID": "CatchID",
        "Length": "Length",
        "Fnode": "Fnode",
        "Tnode": "Tnode",
        "CatchID2": "CatchID2",
        "OrigOID": "OrigOID",
        "Shape_Length": "Shape_Length"
    },
    "geometryType": "esriGeometryPolyline",
    "spatialReference": {
        "wkid": 102113
    },
    "fields": [
        {
            "name": "OBJECTID",
            "type": "esriFieldTypeOID",
            "alias": "OBJECTID"
        },
        {
            "name": "DissolveID",
            "type": "esriFieldTypeInteger",
            "alias": "DissolveID"
        },
        {
            "name": "CatchID",
            "type": "esriFieldTypeDouble",
            "alias": "CatchID"
        },
        {
            "name": "Length",
            "type": "esriFieldTypeDouble",
            "alias": "Length"
        },
        {
            "name": "Fnode",
            "type": "esriFieldTypeInteger",
            "alias": "Fnode"
        },
        {
            "name": "Tnode",
            "type": "esriFieldTypeInteger",
            "alias": "Tnode"
        },
        {
            "name": "CatchID2",
            "type": "esriFieldTypeInteger",
            "alias": "CatchID2"
        },
        {
            "name": "OrigOID",
            "type": "esriFieldTypeInteger",
            "alias": "OrigOID"
        },
        {
            "name": "Shape_Length",
            "type": "esriFieldTypeDouble",
            "alias": "Shape_Length"
        }
    ],
    "features": [
        {
            "attributes": {
                "OBJECTID": 1,
                "DissolveID": 27595,
                "CatchID": 27595,
                "Length": 2781.07146565,
                "Fnode": 26156,
                "Tnode": 26157,
                "CatchID2": 27595,
                "OrigOID": 26064,
                "Shape_Length": 3396.52480192211
            },
            "geometry": {
                "paths": [
                    [
                        [
                            -12199216.57915191,
                            4173685.6358526573
                        ],
                        [
                            -12201719.387905361,
                            4170671.685971223
                        ],
// (snip)
                        [
                            -12201768.97746418,
                            4170677.9241078123
                        ]
                    ]
                ]
            }
        }
    ]
}


texts: Changing shapefile longitude from -180 - 180 to 0 - 360?

texts:
I have a shapefile of world countries, with longitude between -180° and 180°.
How can I change the shapefile (eg. with ogr2ogr) for having longitude between 0 and 360°.

Following OSGEO/Trac explanation about GenParms, I tried this without success (the result is still in -180 - 180:

  ogr2ogr -wrapdateline -t_srs '+proj=latlong +datum=WGS84
  +lon_wrap=-180 +over' /data/tmp/new_algo/world_new.shp /data/tmp/new_algo/world.shp

I need to have the new longitude defined as:

  if Lon < 0:
      new_Lon = 360 + Lon
  else:
      new_Lon = Lon

Of course, the process must take care of the date line.


texts: Multiple reclassification using RasterCalc in QGIS

texts:
Is there any way to perform multiple reclassification using RasterCalc?

For example reclassifying value 1 to 100 in map Map can be written as:

eq([Map]@1, 1, 100)

But if there are multiple values to reclassify in one operation, how can I do it?
For example reclassifying value 1 to 100 AND 2 to 200 in map Map.

eq([Map]@1, 1, 100), eq([Map]@1, 2, 200) //<-This does not work!!

Please give me a advice.


texts: Is it possible to reclassify rasters manually?

texts:
Is it possible to change raster cell values individually or as a group?

Currently, I am digitizing around cells I want to change and then doing a polygon to raster conversion and then a merge to new raster.

I'm using ArcGIS 10.3 for Desktop.


texts:
I have the following query:

SELECT row_to_json(fc)
FROM (
    SELECT ST_AsGeoJSON(lg.geometry)::json As geometry,
            (
                select row_to_json(t)
                from (
                    select name, type, admin_level
                ) t
            )
        As properties

    FROM osm_admin As lg
)

As fc;

This returns all data in table as GeoJSON.

How can I additionally restrict the results to only within the bounding box? (As the second image below.)

I need to get full GeoJSON for a tile. The PostGIS data is probably in 3857, as it's OSM imported as pbf with imposm3.

Have tried this, as per Mark C's suggestion:

SELECT row_to_json(fc)
    FROM (

        SELECT ST_AsGeoJSON(lg.geometry)::json As geometry,

        (SELECT row_to_json(t)
            FROM (SELECT name, type, admin_level)
            WHERE the_geom && ST_SetSRID(ST_MakeBox2D(ST_Point(-73.9980, 40.726), ST_Point(-73.995, 40.723)), 4326)
            As t
        )

        As properties

        FROM osm_admin As lg
    )
As fc;

But I get the following error:

ERROR:  subquery in FROM must have an alias
LINE 6: FROM (SELECT name, type, admin_level)


texts:
So i recently decided to switch my project from Openlayers 3 to Leaflet, mainly because i needed the layer switcher control.

I followed the example on the leaflet tutorial, this one to be precise and i was very happy with the result, however, i really need to "customize" the layer switcher panel, this one i mean:

I want to create my own custom bar and use this to switch the layers, same functions, just different appearance.

I mean, just create a div with radiobuttons/lists on HTML and use this one to switch the layers.

How would i go about doing that? Is there a plugin or something else that can help me?


texts: How to extract GPS information from Garmin Cam 20 video files

texts:
I just got a Garmin Dash Cam 20 and recorded some videos with it, and also confirmed that Garmin's own player can show the GPS information on a map.

I can't seem to find any technical information on where the GPS information is stored (assuming within the AVI file at this point) in contrast to some other Garmin action cam I have that stores the GPS separately in a GPX file. I need to integrate with a system, so manually using Garmin's player to extract it is not an option.

Does anyone have some pointers on where the GPS information might be stored?

There is also audio on the video, so it do not look like its encoded in the audio stream.


texts: Count number of different points within all polygons?

texts:
I have 2 PostGIS tables in PostgreSQL: polytable and pointtable. The point table has 1000s of points and a class column with 6 different unique classes (integer values from 1 to 6). The polygon table has 1000s of polys and a name column. Some names appear multiple times.
For every single polygon, I want to count how many points (per point class) are in each. Although the polygons contain a name column (we'll call polyname here), I don't need these grouped at all... The current query I use only works when using GROUP BY polyname on the polygon table. Therefore I can't treat each polygon individually:
SELECT
  polyname,
  count(pointclass) FILTER (WHERE pointclass = 1) AS pointclass1,
  count(pointclass) FILTER (WHERE pointclass = 2) AS pointclass2,
  count(pointclass) FILTER (WHERE pointclass = 3) AS pointclass3,
  count(pointclass) FILTER (WHERE pointclass = 4) AS pointclass4,
  count(pointclass) FILTER (WHERE pointclass = 5) AS pointclass5,
  count(pointclass) FILTER (WHERE pointclass = 6) AS pointclass6
FROM
    (select poly.polyname, pnt.pointcolumn
    from schema.polytable poly
    left join schema.pointtable pnt
    on st_intersects(pnt.geom, poly.geom)) sub
GROUP BY polyname

How can I remove the grouping in the query above? Simply removing the last line throws an error (column "sub.polyname" must appear in the GROUP BY clause or be used in an aggregate function)
This is my desired output:
polyid   polyname   polygeom   pointclass1   pointclass2   pointclass3...
1        red        0101000..  55            10            5
2        blue       0101000..  3             25            46
3        green      0101000..  100           39            11
4        yellow     0101000..  0             0             3
5        red        0101000..  0             11            37
6        green      0101000..  58            99            200
7        green      0101000..  0             0             0

Note above, how I still want polygon rows returned which have no points within. So a LEFT JOIN of point count to the polygon, per point class. Similar to the query below, but with 6 point count columns for each point class:
SELECT poly.id, poly.name, count(point.shape) AS pointcount
FROM polytable poly
LEFT JOIN pointtable point ON st_intersects(poly.geom, point.grom)
GROUP BY poly.id

I'd prefer not to take the longer route of adding 6 empty columns onto the polygon table for each point class, then running an UPDATE on these columns to populate them, using different WHERE pointclass = 1/2/3... etc... There must be a better way?


texts: Is it possible to show our custom message in QGIS Status bar using python?

texts:
Is there any way to show our custom message in Qgis Status bar using python? Just like in arcgis IApplication.statusbar.message(0) = "Please wait..."

like that is there any option to show progressbar in Qgis like IApplication.progressbar.show()


texts:
I have been trying and failing for some time to convert the radiance to refectance by following step 4 in this tutorial
I have tried everything I can think of but my knowledge is limited. Would anybody be kind enough to please point out how I can fix/avoid the error messages shown in the image here...

I have copied the messages too:

Executing: RasterCalculator rastercalc3_b1_ref=(3.141592654*["rastercalc3_b1"]*Square(1.00679))/(1997.0*Sin(47.40016008*3.141592654/180.0)) \\fsr-srv-06.qm.ds.qmul.ac.uk\Share\RedirectedFolders\fa14012\Desktop\PES\1_data\arc\b1_ref
Start Time: Fri Nov 04 19:58:25 2016
rastercalc3_b1_ref=(3.141592654*[Raster(r"rastercalc3_b1")]*Square(1.00679))/(1997.0*Sin(47.40016008*3.141592654/180.0))
ERROR 000539: Error running expression: rcexec()
Traceback (most recent call last):
  File "<expression>", line 1, in <module>
  File "<string>", line 5, in rcexec
TypeError: can't multiply sequence by non-int of type 'float'

Failed to execute (RasterCalculator).
Failed at Fri Nov 04 19:58:28 2016 (Elapsed Time: 2.38 seconds)


texts: Download quota for Open California dataset?

texts:
I'm trying to download all the available Planet imagery for a site in California by cropping a subarea. I can see there's 112 assets of type REOrthoTile or PSOrthoTile for this site using the stats endpoint. I follow the procedure in the docs about activating & downloading assets, yet I can only download 40 of the 112 images. I get the following error:

{"field": {}, "general": [{"message": "Download quota has been exceeded."}]}

The documentation talks about rate limits, but I can't find any information on the download quota:
https://www.planet.com/docs/reference/data-api/api-mechanics/#rate-limits

The question is -- how to find out whether/what download quotas apply? Is that something one can get programmatically by querying the API (so to be used in an automated client that is mindful of the quota)?

More info on the asset I'm trying to download:

AOI:
{"coordinates": [[[-117.7070372256684, 34.08932328386836], [-117.69617937433158, 34.08932328386836], [-117.69617937433158, 34.08033071613165], [-117.7070372256684, 34.08033071613165], [-117.7070372256684, 34.08932328386836]]], "type": "Polygon"}

item ID: 20151205_192708_1154812_RapidEye-2

item type: REOrthoTile

asset type: analytic


texts:
I've seen similar questions here and I reckon this is a common problem ever since Mac OS High Sierra, but none of the solutions proposed here and elsewhere seemed to work for me. QGIS, after installed, simply won't open. The app logo shows up in the tool bar below when I click, but that's it.

I've tried opening with a right click, reinstalling multiple times (and versions), and changing the security preferences to allow apps from anywhere to be opened, but have found myself with the same results.

Also, I think it might be worth mentioning I had a different MacOS version until a couple hours ago, but when I tried installing QGIS it said I needed a newer version of the MacOS to run it, so I upgraded it to the most recent one (13 Catalina).


texts: I wanted to perform atmospheric correction using MODIS surface reflectance products for validating my atmospherically corrected images (from Landsat 8). As Landsat 8 only became available in 2013, is it safe to use MODIS data of 2013 and beyond?


texts:
I am adding Data from a File Geodatabase to QGIS. Everything is working fine, but the drawing of a specific Line subtyp (first image).
It looks like, QGIS can not draw the circle properly. An export to shape (second image) and adding to Qgis shows how the feature should be drawn.

Anybody seen something similar and a solution?


texts: Label states in QGIS with county-level shape

texts:
I want to label state names in QGIS but my original shapefile is built in county-level. If I simply label them with State, my map shows multiple labels in each state as my original file (one label for each county).

Is there any way to to that in QGIS?


texts: Cartodb; Open an image and pdf file in a new tab from infowindows

texts:
I would like to add a link to open an image and a PDF file in a new tab from an infowindow. Actually, I can open websites from them but not images or PDF files.

For images, I have uploaded them to Imgur and I have tried with the share and direct link but they don't work.....

For PDFs this doesn't work either. I have my file into dropbox (I only know the download link). Where can I upload it in order to show or download it?

Code for the image:

   <b style="color:blue;font-family:Comic Sans MS ;font-size:150%"

     > Link vacio

The PDF file dropbox link:
https://www.dropbox.com/s/pvfymix6sw9ydut/fce_hb_dec08.pdf?dl=0


texts: Openlayers 3 - georeferenced ImageStatic

texts:
I am trying to display a georeferenced image (in 4326) over a Bing map in OL3. The ImageStatic is the following.

var view = new ol.View({
            center: [0, 0],
            zoom: 15,
            projection: 'EPSG:4326'
        });

var imageExtent = [26.545242858515152, 39.083929396585, 26.592433689580588, 39.11195443968136];
var imageProjection = 'EPSG:4326';
var is = new ol.layer.Image({
                    extent: imageExtent,
                    source: new ol.source.ImageStatic({
                        url: 'test/test.tif',
                        imageExtent: imageExtent,
                        projection: imageProjection,
                        imageSize: [8869, 5267]
                    })
                });

The Bing is a usual one:

var Bing = new ol.layer.Tile({
                    source: new ol.source.BingMaps({
                        key: '...',
                        imagerySet: 'AerialWithLabels'})
                    });

and the map:

var map = new ol.Map({
    controls: ol.control.defaults().extend([
        new ol.control.FullScreen(),
        mousePositionControl
    ]),
    interactions: ol.interaction.defaults().extend([
        new ol.interaction.DragRotateAndZoom()
    ]),
    // Use the canvas renderer because it's currently the fastest
    target: 'map',
    view: view
});

The problem is that I cannot see the image. If it is georeferenced correctly shouldn't it be on the right area?


texts: How to get OpenLayers Plugin to render base maps at a large scale i.e. 1:1000 or larger?

texts:
Am I the only one with this problem? You have a nice OpenLayers base map and get ready to digitize vegetation at a very high resolution but as soon as you zoom in, OpenLayers fails to render the base map! It's so frustrating. I've tried using GeopapaTile to create a cache that I can use directly from my hard drive but I can't figure out how to call the zoomed png files (if they are in fact higher resolution, I don't know) that make up the tile cache. I am only able to add the raster overview file that GeopapaTile gives with a tif extension and tfw but when I zoom in, I get pixels, not higher resolution.

I'm using Windows 10 QGIS 2.10.1 OpenLayers Plugin 1.3.6 and both Bing and Google Maps Aerial and Aerial Hybrids. There is another ticket that refers to this problem (http://hub.qgis.org/issues/6822) Question followup "... I have similar issues with Bing aerial, only a bit more severe. In this case, Bing aerial does not display any imagery at all past a certain zoom point." Answer  that says: "The latest issues with Google Maps layers are fixed in the newly released OpenLayers Plugin v1.3.6" (ticket closed)


texts: How can I determine which polygons contain other polygons?

texts:
I am using only free, open source tools. I have two shapefiles: a shapefile of voting precincts and a shapefile of census tracts. I want to know: which voting precincts fit inside which census tracts? What is the mapping vocabulary describing what I am trying to answer?

I am trying to know which polygons from the precinct shapefile are bounded by which polygons from the census shapefile? Is that right? How do I answer that question programmatically? What tools would I use? I have a vague sense that I could do this in PostGIS -- but don't really know. Also, is it possible in QGIS?


texts:
I'm trying to use NOAA's VDatum to verify the NAVD88 -> WGS84 datum conversion in a c++ datum conversion library that I've written.

The answers didn't agree, so I attempted to check a simple ITRF2008 -> ITRF2005 conversion, which should just use a single helmert transformation, and again the results don't agree.

Our helmert transform in question was written using the 7 parameter method in EPSG guidance 7.2 (using the position vector transform notation, the same as IERS), and tested against the GIGS dataset.

Using the (7 parameter) 2008 to 2000 Parameters from IERS Conventions (2010) page 41 yields the following transform*:

X:  1513460.90425574 ->  1513460.90372061
Y: -4463118.57366578 -> -4463118.58057589
Z:  4283412.81104103 ->  4283412.81065335

(* small angle approximation, position vector rotation)

but the VDatum results are several centimeters offset in Z and several millimeters offset in X/Y, which I wouldn't expect if we were using the same transform and parameters:

So my question is, does anyone know whether VDatum uses the IERS specified transformation parameters for ITRFXX to ITRFYY conversions, or even better, what parameters at what epochs are being used?

NOTE: I realize the precision shown here is substantially greater than the accuracy in the parameters themselves, but I'm just trying to verify the math is done correctly, not that the results are accurate to the given number of decimal places. The same transform with the same parameters should match much closer than the second decimal place.


texts: Call a QgsMapToolIdentifyFeature mapTool with a button in a QGIS3 plugin

texts:
In QGIS 3.12, I'm trying to get a feature on the mapCanvas from a plugin trough a pushButton. I used the code from here but it doesn't work in the plugin although it works fine in the QGIS python Editor.

from qgis.gui import QgsMapToolIdentifyFeature

def onFeatureIdentified(feature):
    fid = feature.id()
    print ("feature selected : " + str(fid))

layer = iface.activeLayer()
mc=iface.mapCanvas()
mapTool = QgsMapToolIdentifyFeature(mc)
mapTool.setLayer(layer)
mc.setMapTool(mapTool)
mapTool.featureIdentified.connect(onFeatureIdentified)

I tried to insert it in my plugin:
Connecting the button:
self.dlg.pbAddFeature.clicked.connect(self.addFeature)

def onFeatureIdentified(feature):
    fid = feature.id()
    print ("feature selected : " + str(fid))

def addFeature(self):
    layer = myLayer
    mc=iface.mapCanvas()
    mapTool = QgsMapToolIdentifyFeature(mc)
    mapTool.setLayer(layer)
    mc.setMapTool(mapTool)
    mapTool.featureIdentified.connect(self.onFeatureIdentified)

But nothing happens.
The only way I found to have an action was to extract the maptool as a property in my plugin.
Class MyQGISPlugin
    mapTool = None

then connecting the onFeatureIdentified method in the run method:
self.mapTool = QgsMapToolIdentifyFeature(self.iface.mapCanvas())
self.mapTool.featureIdentified.connect(self.onFeatureIdentified)

then defining the onFeatureIdentified method:
def onFeatureIdentified(feature):
    print(type(feature))

and then defining the press Button method:
def addFeature(self):
    self.mapTool.setLayer(self.myLayer)
    self.iface.mapCanvas().setMapTool(self.mapTool)

When In click the pushButton on my dialog, the cursor changes and I can click in myLayer but it prints me <class 'my_qgis_plugin.my_qgis_plugin.MyQGISPlugin'> instead of a feature. And it prints AttributeError: 'MyQGISPlugin' object has no attribute 'id' when I use the first onFeatureIdentified to print the feature's id.


texts: What is the dimension used in OpenLayers stroke width?

texts:
In OpenLayers, when you are giving a style to layer, say a WFS layer, we can give stroke width, like this:

var layerWFS = new ol.layer.Vector({
    source: sourceWFS,
    title: 'plot_layer',
    name: 'plot_layer',
    style: new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: 'rgba(255, 0, 255, 1)',
            width: 1
        })
    })
});

Here width is 1 , I would like to know in which dimension that we are giving the stroke width, whether:

millimeter (mm) 1 means 1mm
centimeter (cm) 1 means 1cm


texts: Best way to copy GDAL metadata from one JP2 image to another

texts:
I have a separate question going regarding the poor performance of GDAL CreateCopy() on a large (4096x4096) jp2 file, so here's a look at a different angle:

I have a C++ openCV program that takes a JP2 image and makes a copy after doing some manipulation.  The new JP2 image file is missing the GDAL metadata (see example below).  Would it be faster to load the GDALDataSet from the original image, and then call GDAL SetMetaData() for each Metadata item?

If so, does anyone have sample C++ code for iterating through the metadata items in the original GDALDataSet?

Original GDADINFO metadata:

Driver: JP2OpenJPEG/JPEG-2000 driver based on OpenJPEG library
Files: mask.jp2
Size is 4096, 4096
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-97.025756835937500,32.398681640625000)
Pixel Size = (0.000001341104507,-0.000001341104507)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -97.0257568,  32.3986816) ( 97d 1'32.72"W, 32d23'55.25"N)
Lower Left  ( -97.0257568,  32.3931885) ( 97d 1'32.72"W, 32d23'35.48"N)
Upper Right ( -97.0202637,  32.3986816) ( 97d 1'12.95"W, 32d23'55.25"N)
Lower Right ( -97.0202637,  32.3931885) ( 97d 1'12.95"W, 32d23'35.48"N)
Center      ( -97.0230103,  32.3959351) ( 97d 1'22.84"W, 32d23'45.37"N)
Band 1 Block=1024x1024 Type=Byte, ColorInterp=Red
  Overviews: 2048x2048, 1024x1024, 512x512, 256x256
  Overviews: arbitrary
Band 2 Block=1024x1024 Type=Byte, ColorInterp=Green
  Overviews: 2048x2048, 1024x1024, 512x512, 256x256
  Overviews: arbitrary
Band 3 Block=1024x1024 Type=Byte, ColorInterp=Blue
  Overviews: 2048x2048, 1024x1024, 512x512, 256x256
  Overviews: arbitrary


texts: Refreshing QGIS after updating MS Access database?

texts:
Why is QGIS not refreshing after updating a connected database?

It is a MS Access Database which I connected to QGIS using ODBC. The database contains coordinates of some vessels and they are already visible in QGIS but it is only refreshing when I disconnect and then reconnect the db connection.
I am looking for a solution the past few days but I am not possible to find any.


texts: Not able to figure out what is wrong with geolocation or projection in upoaded raster

texts:
I am working on VIIRS SDR dataset. I had different layer of Brightness temperature, latitude and longitude in given .HDF5 dataset. Using reference from http://geoinformaticstutorial.blogspot.in/2014/06/reading-and-map-projecting-amsr-2-data.html I had created different VRT files and warped it to give geolocation. I got raster of  Brightness temperature as shown in uploaded image below. it is in .tif format. I am also attaching VRT files which I had created. https://drive.google.com/open?id=1t3gFJYl0SdysS2zEJxlRJsU0kD0CR0_J (Bright.vrt) , https://drive.google.com/open?id=1hbCyWGjqw4VwIJjk9d_bQKW1ElJ5ibxb (lat.vrt) and https://drive.google.com/open?id=1ni27bXVnQ0ejyHMmF47Hebba-FOFTP3x(lon.vrt)
According to user guide of the dataset, lines appearing in the image should be removed after giving geolocation.

Kindly help me to figure out what can be wrong with it?


texts: Is there a way to cut polygons in QGIS using freehand tool? Unfortunately Freehand Editing tool in QGIS 2.x does not provide such functionality. There is a tool for cutting polygons with standard editing mode . But I need to cut polygons with freehand mode.Is it possible?


texts: How to manipulate select query in mapinfo so that it does not high light the selected data?

texts:
I have select query in a button but I don't want user to see the way I am changing the map coordinate while selecting specific points and adding a layer. How is it possible to do that in a button ?
 Is using select query the correct way or are there other ways?


texts: How do I make this into a working CRS in ArcGIS?

texts:
+proj=tmerc +lat_0=0 +lon_0=-114 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

I have this text string which I successfully used to project data in QGIS. How do I use this to define a CRS in ArcGIS and project my data correctly? Thanks


texts: Finding nearest raster cell value based on vector point?

texts:
I have two layers in ArcGIS: One is a vector point and the other is a raster layer. I want to append the value of the raster data grid to the vector point. The problem is the vector point is away from any existing raster grids with value, so the command "Sample" return no value for those vector point.

So my question is: How can I find the nearest raster cell given a vector point and extract the value from the cell?

Actually there are two things I need to do:

Identify the outliner within a ranger (Usually 1-2 cells away from my existing raster)
Assign a value from the nearest raster to them, if they are in a range

Edited:
I have about 3000 vector points and my problem is this:

I need to append the raster value based on locations to the vector points. I used the tool "Sample" and it works well for most of the points.

However, there is a situation like this

The right point can get the value from the raster ("Sample" works) but the left point cannot due to the alignment issues.

Raster to polygon is not working because I get error message that I cannot turn them into vector because they are out of the domain. Also I have points like this

Those upper points should not get any values because they are too far away from the raster.

I thought of someways to work around:

Step 1. Do the "sample" first

Step 2. Select the null value after the "sample" operations

Step 3. construct buffer zones (within 1 decimal degree) based on the null points

Step 4. use zonal statistics??? or other tools to loop through all the raster cells in the buffer, find the nearest one, extract its value and put it in the vector point.

I am stucked at step 4. I don't know which tool in ArcSDK I can use to achieve this function.

Or else...Do you have better idea to deal with that?

I have 20 raster layers like this and I would like to make it automatically (using model builder and arcobject).


texts:
I am trying to manipulate bad pixel values in a raster by using the Con-Statement via ArcPy. Using the build-in Map Algebra function "Raster Calculator", the results are good. To transfer the work flow to a script, the outcome is different. The pixel values should be between [-1;1].

arcpy.env.workspace = out_folder #path to my folder containing the raster files
list = arcpy.ListFiles('*.tif')
for ras in list:
     outras1 = Con(ras > 1, 1, ras)
     outras2 = Con(outras1 < -1, -1, outras1)
     outras2.save(ras + '_mod.tif')

Can you tell, what I am doing wrong?

The results show only pixel values of 1. Additional information: I have large areas of NoData values in my raster.


texts: Making rhomboid with user defined angle in ArcMap?

texts:
I need to capture some buildings in ArcMap 10.3.

It works fine for perfectly rectangle buildings, but whenever some location has buildings with an angle like a rhomboid then it became very time consuming to first draw a rectangle and then edit its vertices to fit.

In the picture above blue ones are the correct building shape and green one is default rectangle shape. Note the angle difference.

Is there any extension, add-in, by which a user can define an angle first and whenever he makes a rectangle then it automatically edit its vertices to fit that particular angle?


texts:
I have a few services added as layers in a map, that is configured to be used in ArcGIS Collector(20.2.1), and I am facing an issue with the synchronization of my edits.
I checked the REST end of the service and found that the replicas are no more available and have disappeared.
I have isolated one of the sync errors with a couple of services and it seems that the corresponding replicas have disappeared from the list of registered replicas:
{"error":{" code":404, "message": "Unable to synchronize replica.", "details":["Replica with GUID: 0cc27ea0-0391-455f-8432-70e3fca794c9 does not exist on the server."]}}

What could be the reason for the disappeared replica?


texts: Calculating yearly forest gain from Hansen data set

texts:
I want to calculate yearly forest gain in my specified area from the Hansen data set - I used the code given by google earth engine online and replaced the word 'loss' with 'gain' but it won't work saying 'gainyear' does not match any bands. How do I fix this?

This is my code:

var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015');
var gainImage = gfc2014.select(['gain']);
var gainAreaImage = gainImage.multiply(ee.Image.pixelArea());

var gainYear = gfc2014.select(['gainyear']);
var gainByYear = gainAreaImage.addBands(gainYear).reduceRegion({
  reducer: ee.Reducer.sum().group({
    groupField: 1
    }),
  geometry: tapajos,
  scale: 30,
  maxPixels: 1e9
});
print(gainByYear);

var statsFormatted = ee.List(gainByYear.get('groups'))
  .map(function(el) {
    var d = ee.Dictionary(el);
    return [ee.Number(d.get('group')).format("20%02d"), d.get('sum')];
  });
var statsDictionary = ee.Dictionary(statsFormatted.flatten());
print(statsDictionary);

var chart = ui.Chart.array.values({
  array: statsDictionary.values(),
  axis: 0,
  xLabels: statsDictionary.keys()
}).setChartType('ColumnChart')
  .setOptions({
    title: 'Yearly Forest Gain',
    hAxis: {title: 'Year', format: '####'},
    vAxis: {title: 'Area (square meters)'},
    legend: { position: "none" },
    lineWidth: 1,
    pointSize: 3
  });
print(chart);


texts: Atlas with different projections in QGIS

texts: I have created an atlas -- i.e. a map with a layer defining some regions of interest -- in QGIS. I can render maps of each region using the "atlas" feature of the print composer, but I want to show different regions with different projections. Is that possible?


texts:
I have some XYZ bathymetry data in an irregular, non north-up pattern -- i.e. not gridded. It is saved as a CSV file.

443564.526,6636541.778,-105.47
443564.526,6636538.778,-105.44
443567.526,6636544.778,-105.54

I can import the data into QGIS as Delimited and save as a Vector.
What I want to be able to do is use the various Contour, slope analysis, Shaded relief etc tools.
What do I need to do? Most of the answers I have found need either gridded data or externally generated raster files.

I have found this, but I don't know how I go about using it.
https://www.gdal.org/gdal_grid.html


texts: How to calculate the difference between two rasters in QGIS?

texts:
I have two similar rasters of the Mount Vesuvius and surroundings: one spans 20 kilometers, and the other one 50 kilometers. The two rasters are overlapping:

What I would like to do is to create a new raster file which will show the difference in elevation between these two rasters. I tried using the Raster calculator with the following expression:

"vesuvius_radius_20KM@1" - "vesuvius_radius_50KM@1"

But for some reason it does not work. After I press OK, all I get is some sort of black raster file which does not have any kind of elevation data in it, nor it has any number of pixels:

What am I doing wrong, and how can I create a raster file which will show the difference between the upper mentioned two raster files?

Here are those two raster files:

https://www.dropbox.com/s/98hcftjsnmrqs8p/vesuvius_radius_20KM.tif?dl=0
https://www.dropbox.com/s/a68921e6tpszt0f/vesuvius_radius_50KM.tif?dl=0


texts: Wrong or missing parameter QgsProcessingParameterFile QGIS 3.4 Python plugin

texts:
I'm making a script that can be called by the processing.execAlgorithmDialog function. For one of the parameters, I want to use the QgsProcessingParameterFile class for selecting a .csv file.

I tried the following:

class T4_Visualiseer_Networks_Resultaat(QgsProcessingAlgorithm):

    PrmOutputFile = 'Networks result'

    self.addParameter(
        QgsProcessingParameterFile(
            name=self.PrmOutputFile,
            description="Networks resultaat",
            extension="csv(*.csv)"
        )
    )

    def processAlgorithm(self, parameters, context, feedback):
        # Some code to do after clicking ok

When I select a .csv file in a folder and press OK to begin the script, I get the message Wrong or missing parameter value: Networks resultaat (see picture below).

I already tried the following:

Replace / by \
Add behavior=0 to make sure a file is selected (this was already the case because the default value for behavior is 0)
Rename the .csv file

None of the above actions fixed the problem of getting this message and I can't find any other solutions on the internet.

Do you have any idea what the problem can be or how I can fix it?


texts: Package Layer (Data Management) tool reporting error 000814?

texts:
I want to package layers using the Package Layer (Data Management) tool but I always get error 000814 indicates invalid output file type (I try to add the .lpk extension also the same error).

The following pictures show the processing of package layer and the error message.

I search internet and find it is in ArcGIS 10.2 Issues Addressed List (My production version of ArcGIS is version 10.1).

I want to confirm is it a bug or is there any means to compress layers (in enterprise geodatabase) for both 10.1 and 10.2?

In ArcMap “Table Of Contents”, I select the layers to be packged, right click and choose “Create Layer Package” ->save package to file -> navigate the a specified location -> click "save".

After that I use windows browser and arccatalog to access the *.lpk, but I cannot find the package file in the specified location.


texts: Unable to Use Defined Variable for Raster Calculation in Python

texts:
I am trying to open a text file that contains one line, which is a number. I am then trying to define this as a variable to use in my raster calculation/map algebra. I am currently running on Windows 7.

Command Prompt Error

RuntimeError: ERROR 000732: Input Raster: Dataset slope1 does not exist or is not supported

input2 File

-.0000462

Script Example

slope1 = open('C:\Users\Apples\Documents\ArcGIS\Blueberry_Cobbler\input2', 'r').read()

outraster = ((slope1 * Power(inraster ,3)) + (slope1 * Power(inraster ,2)) + (0.0056 * inraster)  + 65.864)


texts:
I use QGIS 2.18.14
I import csv file with hits on a certain area. I then create a simple indicative heatmap (properties > heatmap).

Can I add contour lines to this heatmap and if yes, how?

I also put a vector grid on these hits > count points in polygon and then via properties > graduated based number of hits per square.

Can I add contour lines to this square based heat map and if yes, how?


texts: How to create zip code areas from point data?

texts: I have a point dataset which has zipcodes for a certain region. This zipcodes have been marked on the basis of lat long collected. So each zipcode is repeated as it would be useful to show the region covered by one zipcode. Now since there are points I would like to know if there is any tool in QGIS or ArcGIS which can make polygons for each zipcode.


texts: Error importing Shapefile to PostGIS using Importer GUI

texts:
I'm trying to import a shapefile into a PostGIS database using the graphical importer. I just installed POSTGRESQL, and PostGIS. The GUI looks like this:

After testing that my database worked, I chose the shapefile that I wanted to upload, and then received a message that it had failed to import. I tried tweaking a few paramaters, but I'm not sure what I'm doing.

Are there any good tutorials or instructions about how to do this? I'm using Windows 7 64 bit.

Just adding a little more of my error message.

  Shapefile import failed. Connection:
  host=localhost port=5432 user=postgres
  password='******'  Destination:
  public.Area Source File:
  D:\files\Area
  Shapefile type: Polygon Postgis type:
  MULTIPOLYGON[2] Failed SQL begins:
  "SET CLIENT_ENCODING TO UTF8; SET
  STANDARD_CONFORMING_STRINGS TO ON;
  BEGIN; CREATE TABLE "public"."Area"
  (gid serial PRIMARY KEY, "fid_1" int4,
  "area" numeric, "dtm" float8, "dsm"
  float8, "hgt" float8, "nat_area"
  numeric, "nat_peri" numeric, "nat_vol"
  num" Failed in pgui_exec(): ERROR:
  type "geography" does not exist LINE
  14: "the_geom"
  geography(MULTIPOLYGON,4326));
                      ^

  Shapefile import failed.

I've tried running @MerseyViking's suggestion, and then importing again. In the Geometry Column: box I chose MULTIPOLYGON and then tried to import the projected shapefile. The following errors appeared:

Failed SQL begins: "CREATE INDEX "NeighborCheck_MULTIPOLYGON_gist" ON "public"."NeighborCheck" using gist ("MULTIPOLYGON" gist_geometry_ops);
COMMIT;
"
Failed in pgui_exec(): ERROR:  current transaction is aborted, commands ignored until end of transaction block

Shapefile import failed.

Then I tried ticking the Load into GEOGRAPHY column box in Options, and received the following error:

Failed SQL begins: "CREATE INDEX "NeighborCheck_MULTIPOLYGON_gist" ON "public"."NeighborCheck" using gist ("MULTIPOLYGON" gist_geography_ops);
COMMIT;
"
Failed in pgui_exec(): ERROR:  current transaction is aborted, commands ignored until end of transaction block

Shapefile import failed.

Both error messages are identical and I don't quite understand them.


texts: ArcGIS Collector Sync Replica Internal Error?

texts:
I have downloaded maps on arc collector for offline use but whenever I try to sync, it gives me the message "Unable to synchronize replica. Sync Replica Internal Error". To note, many users have been using the same map from different devices.

To experiment, I created a new map with the same feature layer as the map mentioned above, downloaded it to a device, created a point, and had no issues with syncing.

Long story short, I am assuming that my issue is that I have several users adding data to the same map. Everyone downloaded the map to their device since we are using this for emergency response after a hurricane and were not sure if we would have cell service after the storm. Does anyone know of a solution to this problem?

If necessary, I can transcribe the data from each device individually, but that is not ideal.


texts: Drawing Rectangle with INewRectangleFeedback?

texts:
I´m developing in ArcGIS Engine and my question is how to draw a rectangle with INewRectangleFeedback interface. Here is my code but it does not create the rectangle. If someone could help me it would greatly appreciate it.

Public Sub mousedownFBRectanglePoly(ByVal px As Long, ByVal py As Long)

    feedbackStartPoint = New ESRI.ArcGIS.Geometry.Point
    feedbackMovePoint = New ESRI.ArcGIS.Geometry.Point
    feedbackScreenDisplay = Principal.m_mapControl.ActiveView.ScreenDisplay
    feedbackFillSymbol.Color.CMYK = RGB(255, 0, 0)

    'Initialize a new RectangleFB object
    NewRectangleFeedback = New NewRectangleFeedback
    NewRectangleFeedback.Display = feedbackScreenDisplay
    NewRectangleFeedback.Symbol = feedbackFillSymbol 'QI

    'pass the start point from the mouse position, transforming it to map coordinates.
    feedbackStartPoint = feedbackScreenDisplay.DisplayTransformation.ToMapPoint(px, py) 'QI
    NewRectangleFeedback.Start(feedbackStartPoint)
End Sub

Public Sub mousemoveFBRectanglePolygon(ByVal px As Long, ByVal py As Long)
    '    'pass X and Y to feedbackMovePoint to transfer to NewEnvelopeFeedback
    feedbackMovePoint = feedbackScreenDisplay.DisplayTransformation.ToMapPoint(px, py) 'QI
    NewRectangleFeedback.SetPonint(feedbackMovePoint)
End Sub

Public Sub DoubleclickFBRectanglePolygon(ByVal px As Long, ByVal py As Long)
    'feedbackStopPoint = feedbackScreenDisplay.DisplayTransformation.ToMapPoint(px, py)
    'when mouse comes up, end the new envelope and pass it to feedbackEnvelope.
    feedbackRectangle = NewRectangleFeedback.Stop(feedbackStopPoint) 'QI
    'pass the new rectangle element, the geometry defined by our feedback object
    feedbackElement = New RectangleElement
    feedbackElement.Geometry = feedbackRectangle
    'make sure the element is activated in the current view
    feedbackElement.Activate(feedbackScreenDisplay)
    'now add the newly created element to the ActiveView with default symbology.
    Principal.m_mapControl.ActiveView.GraphicsContainer.AddElement(feedbackElement, 0)
    'and refresh the view so we can see the changes.
    Principal.m_mapControl.ActiveView.Refresh()
End Sub

And also I have the problem that when I call my methods in the mouse events it does not enter the method; it skips even if all the parameters we pass are good.


texts: Adding feature layers in Esri's WebAppBuilder without using Add Data widget

texts:
I've added a feature layer to the "my content" section but the only way I can seem to add it to the WebAppBuilder is through the AddData widget.  When I try adding it through the "Attribute -> Extra Data Source" then it inds the feature layer, but nothing displays on the map.  Then, in a custom widget I'm trying to create, I tried to use the Javascript API to add the feature layer to the map object.  It doesn't throw any errors, it just doesn't display the layer.

Any ideas? I'm not sure why the AddData widget would be different from the other methods, I tried looking through the code but didn't find anything that jumped out at me.  Maybe I'm just missing something small?


texts: Adding new WMS Data Source to QGIS using Python

texts:
When using PyQgis for QGIS 3.4 is it possible to add a new WMS/WMTS Server to the Data Source Manager?

If yes and once created, is it then possible to add a new WMS layer to the map using the new connection within a Python plugin?

I have been trying to add a layer using authentication but my server needs to pre-authenticate so that is not working and in any case the only way it works using the normal GUI is to select the two ignore GetCapabilities options in the Data Source Manager anyway so it is unlikely it would work via QgsRasterLayer anyway.


texts: Differences between OGC and ISO WKT/WKB specification?

texts: Are the OGC-Simple feature access and the ISO (ISO/IEC 13249-3:2011 "Information technology - Database languages - SQL multimedia and application packages - Part 3: Spatial") significantly different (different enought to justify purchasing the  ISO/IEC 13249-3:2011 PDF)?


texts: Openlayers 3 + Use React to display vector features

texts:
I am using openlayers in a React based JS app. I have a vector layer that has different colors according to a field given from the server.

Is it possible to use react for drawing the icons on the map? or should I use simple openlayers object with simple plain javascript?

Regards,
Ido


texts: Running database operations against a stored GIS geometry

texts:
I'm using Quantum GIS and Postgres with the GIS extensions to store shapefiles that I'd like to run point in geometry and point in union geometries against. While it's easy to save and load the vector layers to the database, I'm having problems figuring out how I can run operations against them. How do I use the database to SELECT WHERE lat/long values in another table are contained with a saved polygon geometry? How do I union a couple of geometries and use that in the WHERE clause? I'm having a lot of trouble figuring out what the schema of the saved polygons actually is.

Edit: The shapes have a number of holes in them, if it matters for performing the selects and the unions.


texts: Openlayers WFS and Mapfile: Possible to specify multiple online resources and generic URL/variables in OpenLayers file?

texts:
I have a working mapfile with an associated working OpenLayers file.

The mapfile has a wfs_onlineresource like so:

"http://one_of_my_servers:8190/cgi-bin/mapserv.exe?"

Then in the OpenLayers js/html file you have the wfs layer like so:

var wfs = new OpenLayers.Layer.Vector("Airports", {
            styleMap: new OpenLayers.StyleMap({
                    externalGraphic: "red_pin.png",
                    graphicWidth: 18,
                    graphicHeight: 18
            }),
            wrapDateLine: true,
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.WFS({
                srsName: "EPSG:3857",
                version: "1.0.0",
                url: "http://one_of_my_servers:8190/cgi-bin/mapserv.exe?map=/tmp/navhome/navcore/airports_viewer_svn/airports_viewer.map&SERVICE=WFS&VERSION=1.1.0&srsName=EPSG:3857",
                featureType: "Airports"
            })
    });

Like I said this works fine. Now the issue:

we have quite a few servers here. Say you want to take these two above files and go to server 177. Is there a workaround so you don't have to go into the mapfile and specify/change the wfs_onlineresource and go into the OpenLayers file and change the url: in the wfs layer?

Basically I'm trying to make a 'universal' edition so matter what the server, it will work. Not really sure if this is possible after reading the WFS Mapfile Documentation:

See section 12.3.3 of the WFS 1.0.0 specification for the whole story about the online resource URL. Basically, what you need is a complete HTTP URL including the http:// prefix, hostname, script name, potentially a “map=” parameter, and and terminated by ”?” or “&”.

This may be a longshot, but any thoughts on the matter are appreciated.


texts: Cursor gives "missing field" error, but doesn't say which field is missing

texts:
I am using a search cursor as follows:

fields = ["List", "of", "fields"]
FeatureClass = "path to a feature class, shapefile, xls, or dwg"
with arcpy.da.SearchCursor(FeatureClass, fields) as cursor:
    for row in cursor:
        stuff happens

and I am getting an error when searching an xls file:

Traceback (most recent call last):

File "C:\Workspace\Check.py", line 3, in <module>
    for row in cursor:
RuntimeError: No value given for one or more required parameters.

Failed to execute (Check).

And I get a similar error when searching a .shp file or feature class:

a column was specified that does not exist

I know from experience it is trying to tell me that a field in my file is missing, so my question is is there a way in the code to get it to tell me WHICH field is missing?  I get this error when using a search cursor on both .xls files and .dwg files, and I get a similar missing field error on .shp and feature classes sometimes.  I would like a general way for a search cursor missing a field to tell me which field I am missing.


texts: Map labels are smaller in exported image than in print composer

texts:
I am trying to export images with legible town name sizes but they keep shrinking to an unreadable size when I export them. I have tried a few solutions based on similar questions but to no avail.

I am using QGIS version 2.18.9 and a number of different map sources (OSM, SIX, Google, Bing).


texts: Work with MOD07 - Retrieved Moisture Profile layer - data in R

texts:
I downloaded a complete month of daily images of MODIS MOD07_L2 product (Atmospheric profile). I've worked with *.hdf files for really long (using both R and tools such as HEG), nevertheless, I am now stuck with MOD07 product. I need to work with Dew-point Temperature at different pressure levels, but I don't understand how the data is actually saved in the *.hdf layer (which is Retrieved_Moisture_Profile). I tried using HEG software with all possible options (for output file types), but none of them worked. For every iteration I tried in HEG I got the "Could not find TimeofDay attribute, CalendarDate, Sensorname, pointingangle (and many more) attribute in input hdf file". So I believe HEG is not capable to work with this type of data. The result for every image was a black square with just one pixel value.
I opened the *.hdf file directly in QGIS (version 3.10.4-A Coruña) and, more specifically, I opened the Retrieved_Moisture_Profile layer. When I studied the histogram of the latter image it is quite clear that there are 20 different Dew-point temperature profiles (each of one of them is for an specific pressure level), as you can see in the following figure.

I was thinking that maybe I can just save the *.hdf to *.nc, but I tried with rhdf5 and the latter package is depricated for R version 3.6.3. I also tried opening the *.hdf file with gdalUtils, which actually works, but I couldn't find how to 'divide' each layer of the raster or even save them separately.
So, what I need is the following (in R).

I need to open the MOD07 *.hdf file
I need to open the Retrieved_Moisture_Profile layer within the *.hdf file
I need to save the dew-point temperature profile for each pressure level i.e. for each 'band' (from 1 to 20), as a GeoTIFF file.

FYI: I'm running R version 3.6.3 in Windows 10 64 bit machine.
Thanks in advance for your help!
Sincerely, Jorge.


texts: Using Field Calculator to Combine Month, Day, Year into Date

texts:
I merged 2 layers in QGIS and am trying to clean up the attribute table using the Field Calculator. One layer has 1 column called Date in the format YYYY-MM-DD. The other layer has 3 separate columns with real numbers for Month, Day, and Year. I'd like them all to be in the YYYY-MM-DD format, preferably using the Field Calculator. I'm not experienced with Python and don't know where to begin with the expression builder.
So here's an example of each type of date format within the same layer.
Month, Day, Year all in separate columns

And this is how I want all of them to look (YYYY-MM-DD)


texts: Is it possible to have the attributes fields already filled from the previous object?

texts: I need to create a lot of points with lots of attributes, but most of the fields are the same, so I don't want to edit all the fields all the time. Is there a "remember previous fields" function available?


texts: Looking for documentation of OGR SetSpatialFilter

texts:
I need documentation on how to apply a spatial filter (using OGR's SetSpatialFilter).

Setting the feature layer by which another layer should be filtered seems rather straightforward. But, how do you set the type of spatial filtering ('intersect', 'contain' etc.). And where do you find a list of optional types?

In the case at hand, I would like to filter out features from one layer (let's call it 'TilePolys')that are touched ('intersects with' (?)) by features from another layer (let's call it 'Buildings).


texts: Print boundaries of one layer over another

texts:
i am new to QGIS. I have a base map from google, openstreet map over which i want to show the boundaries of another layer. How can i do that?

While trying to edit, is it possible that QGIS allows to lay emphasize on the boundary of one layer over the base map of the google maps?


texts:
I am switching from Mapinfo Pro to QGIS and have a question. In Mapinfo, I have coded a theme selector to turn on/off themes (workspaces) in a basemap, as shown below. Has anyone made a similar plugin for QGIS, that turn on/off QGIS-projects?

Best Regards
Jan Hvingel


texts: Why is my DEM raster (.tif) not converting to TIN (Gives "Error 999999" instead)

texts:
I am working with a DEM raster file (.tif) which I am trying to convert into a TIN. However it won't let me do this and instead I get the following "ERROR 999999: Something unexpected caused the tool to fail. Please refer to the error help for potential solutions, or contact ESRI technological support".

I checked out the error description before but I got lost. What can I do in this case? The file I work with is an 0.5m * 0.5m raster of the area its measured heigh. Furthermore, the file is a bit under 500mb.


texts: Problems with "Filter" and "Select with Expression" statements in QGIS

texts:
I am using QGIS 3.10 to make simple queries using the "Filter" and "Select by Expression" tools. I have a layer with cities/towns in New York State with population attributes, e.g. "POP1990" is population in 1990, "POP2010" is population in 2010, etc.
Certain queries seem to work in one tool but not the other. For example, I want to select all towns/cities in New York State where the population increased by more than 5% from 1990 to 2010. Oddly (to me) when I make all references to an attribute/field name lowercase, the "Select by Expression" tool works, but when I spell the attribute/field name with the correct case and put it in double quotes it does not work.
Below are the expressions I write, the tool I use, and the outcome:
Tool: "Select by Expression"
Expression: ( "POP1990" * 1.05 ) < "POP2010"
Outcome: doesn't work - selects all features
Tool: "Filter"
Expression: ( "POP1990" * 1.05 ) < "POP2010"
Outcome: doesn't work - all features remain
Tool: "Select by Expression"
Expression: ( pop1990 * 1.05 ) < pop2010
Outcome: WORKS CORRECTLY - 469 out of 994 features selected
Tool: "Filter"
Expression: ( pop1990 * 1.05 ) < pop2010
Outcome: doesn't work - all features remain
I could give other examples of different expressions working or not working.
Data are downloadable below as a link, if you would like to test it yourself. (This is my first Stack Exchange post, by the way.)
Here is the New York State cities/towns shapefile layer for download.


texts: Why would a map lost the coordinate system when I use it in ArcGis?

texts: I have a file in AutoCad with ETRS89 coordinates, and when I use the same map in ArcGIS it loses coordinates. With other files do not happen.


texts: Calculating weighted values using function in Function Editor

texts:
After my previous question Weighted values of field in QGIS, I successfully adapted the following formula in Expression dialogue
"Value" / (sum("Value") / (count(CASE WHEN "Value" >= 1 THEN +1 END)))

In QGIS in the Attribute Table, there are many fields with values for which I count using this formula.
However, I don't want to repeat all long action and therefore I need to create a function in the Function editor. How can I implement this?


texts: send SPECS to mapFish service using POST method

texts:
I'm using mapFish print-servlet-1.1 and I generate the PDF from the default interface (index.html) from tomcat.
The print service is a standalone service deploied on tomcat 7.

If the form method is GET, everything is wotking fine.
If I want to use POST (to send more SPECS) I get this result:
Error while generating PDF:
Unknown method: /print.pdf


texts: Steiner tree problem: Optimal routes for Fibre planning on road network

texts:
I am trying to solve a large optimization problem where I have a set of nodes on a road network that need to be connected using fibre using the shortest possible amount of fibre.

The problem is a classical Steiner tree problem. Grass has a function called v.net.steiner but my road network consists of 20,000 edges connecting intersections and the number of Steiner nodes (points to be connected by fibre) is about 1500. I tried this on a 32-bit Linux OS with 4 GB of memory and I got an out of memory error.

Any suggestions on how to solve this with a desktop PC - or do I need to send the problem to a super-computer

I've attached an image of the problem - don't worry I will trim the dangle-roads but that only shaves off about 1000 road segments


texts: Changing ArcSDE binary formatted data to MSSQL Geography or Geometry data types?

texts:
We are currently using ArcSDE 9.3 to store our Geodata. In a project which is currently being implemented, the software being used can access most data formats including the ArcSDE binary format and/or the MSSQL Geometry or Geography format. Our current problem is that the new software cannot implement complicated nested sql queries or use certain functions on the SDE binary format with the ArcSDE 9.3 API. According to the software firm, these features apparently only work with the MSSQL geometry or geography format.

Based on the context above, I am looking for some advice regarding the following:

How can I change all my data which is currently stored as SDE Binary into the MSSQL Geometry or Geography data format?
Are MSSQL Geometry and Geography formats fully compatible with ArcSDE and ArcGIS?
What disadvantages do Geometry and Geography formats have in comparison to SDE binary?


texts:
I have a point data set that is recorded close to Ottawa in Canada. I am going to use it in some training and for the training purposes, I would like to move the data set to a place a bit east of Cape Town in South Africa, i.e I want to move the data approximately 90 degrees eastwards and 80 degrees southwards.

Is there any way to do this in QGIS (I don't mind running a Python script in the console) or do I have to write a Python script to parse the file I initially imported into QGIS?

The data are presently used in QGIS as WGS84 UTM 18N, but I have no problem reprojecting it to WGS84 lat/lon if that makes things easier.

(I know moving the dataset to another latitude will change the shape of it, it may be that I need to be a bit more careful with the CRSs like - moving it to the southern hemisphere in the same UTM zone, then move it eastwards by manipulating UTM zones and adjusting the easting values, but for a first attempt I am pretty sure lat=lat-80, lon=lon+90 would be good enough for what I need it for)


texts: GeoTIFF rendered in OpenLayers using GeoServer: Tiles have different brightness and sharp edges between them

texts:
I have a script that download and process satellite images using pyhon/gdal and then upload them to GeoServer.

The basic process is:

gdalbuildvrt tmp.vrt file_1.tif file_2.tif ... file_n.tif
gdal_translate -of GTiff -co TILED=YES tmp.vrt outfile.tif
gdaladdo -r gauss outfile.tif 2 4 8 16 32

#then upload to geoserver

outfile.tif looks like this in QGIS:

When GeoServer feeds image tiles to my OpenLayers application I am to understand that it applies a contrast enhancement. This is initially a good thing since my images tend to be a little dark.

The problem is that the effect seems to be applied locally on each tile. This is how the same file looks like in OpenLayers:

I have tried to modify the used style file (SLD) in GeoServer. Changing the ContrastEnhancement property from Normalize to Histogram makes the effect even bigger.

How do I remove these "boxes" in the image?

Does someone have any ideas for what I could try?


texts:
I can't seem to figure out or find a work around for the symbology labels after a reclassify when trying to automate maps. I run the following code and it is almost doing exactly what i want it to do.  After I change a data source path I run the following 2 lines of code.  The layer is set to have the colors i want and the numbers match what should be there. The only issue is the display of the labels for the classification.

                arcpy.mapping.UpdateLayer(df, lyr, updateLayer, True)
                lyr.symbology.reclassify()

The map is set to have 1 decimal place as is the layer file but after the reclass it is set to 4 significant digits.

Does anyone have any thoughts?

This is the one thing i need to fully automate the maps i am working on.  It seems crazy that after doing a change path and a reclass that i have to open the mxd reset the layer properties and then update notes that i have generated based on the new classes with the update ranges.


texts: Split By Attribute for point feature gives error, 'unicode' object has no attribute 'polygon' in ArcMap 10.7

texts:
I have a point feature where many points share the same value as other points on an attribute called "site".  I'd like to group them by that attribute, and then make the groups into a seperate feature for each group.
I'm trying the "Split By Attribute", but keep getting the Error,

  'unicode' object has no attribute 'polygon'

Searching that error only gives my info on Python scripts.  I thought maybe it only works on polygon layers, but I'm getting the same error when I try to spit up a layer of polygons.
Does anyone have any idea what I'm doing wrong?

Below are the parameters I'm using:

Input Table

CWIRS Agreggate/CWIRS sites  -  point feature

Target Workspace

Scratch Geodatabase,  I've also tried a directory.

Split Field

Is just a text value


texts: Correctly setting up CRS systems for the NUTS (Eurostat) geographies and Doogal's data?

texts:
I'm interested in undertaking spatial matching in R with use of the sp package on some of the publicly available data consisting of:

Polygons: Nomenclature of territorial units for statistics
X/Y coordinates: sourced from Doogal's website, with Easting values being assigned to x and Northing to y.

Problem

The challenge seems to be with setting up the correct CRS for both of those data sets.

Illustration

I've initially attempted to set up projections to those discussed here, via:

shpNUTS <- spTransform(
  shpNUTS,
  CRSobj = CRS(
    "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
  )
)
proj4string(pointsDta) <- proj4string(shpNUTS)
res <- over(pointsDta, shpNUTS)

results in an unusable map:

Similarly, the attempt to remove projection via:

proj4string(shpNUTS) <- CRS(as.character(NA))

Does not improve the map.

Matching

The matching is undertaken in a simple manner,

sp::over(pointsDta, sp::geometry(shpNUTS))

but I've also tried:

res2 <- spatialEco::point.in.poly(as(pointsDta, "SpatialPointsDataFrame"), shpNUTS)

The matching does not return overlapping points; this is obviously wrong as x/y coordinates sourced from Doogal'swebsite reflect UK postcodes so all of the postcodes should be matched against a set of UK NUTS geographies.

Remaining code

The shapefile is read via the readOGR function available through rgdal package, Doogal's data is read via readCSV in R. Remaining code reflects basic transformations:

pointsDta <- doogalDta[, c("x", "y")]
pointsDta <- na.omit(pointsDta)
coordinates(pointsDta) <- ~ x + y
# reproject nuts
shpNUTS <- spTransform(
  shpNUTS,
  CRSobj = CRS(
    "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
  )
)

proj4string(pointsDta) <- proj4string(shpNUTS)
res <- over(pointsDta, shpNUTS)

Update

Updates in reply to provided comments.

Shapefiles

Concering the shapefiles, I've sourced this file, that is made available through the Eurostat website.


texts: Calculating Aerosol optical thickness from Landsat 7

texts:
EDIT I have added question down
I am reading one paper and I found method to retrive AOT from landsat 7 image
This should be simplifield equation, because they assumed observation zenit as 0. Than they put on landsat image grid cell 20x20.
So have they calculated standard deviation for each cell and than with this equation they should get AOT also for each cell. Am I correct?
EDIT And at the end if I have AOT for first cell I just put in raster calc [raster1] * AOT1


texts: Choosing color classes for comparing differently distributed spatial data

texts:
For comparison, I want to draw two maps alongside each other. Both maps should use the same color classes/scale. The problem is that the data for the maps is very different: The data for the first map (A) consists of many small spatial units while the data for the second map (B) consists of fewer and larger spatial units but also some outliers with regard to the measured attribute.

Now, how I should determine color classes?
If each color class covers the same range (equal intervals) then the outlier of data set B will skew the color scaler, making relevant differences indistinguishable.

If I use jenks intervals or quantiles then the values of data B will be underrepresented because of their smaller number.

What can I do about the problem?


texts: Setting the visual style of a raster from python console in QGIS

texts:
I am trying to give proper format to a raster (from a GEOTIFF file, previously generated). I need to show it with 3 classes, using discrete classification and equal intervals. This is easy to do from within the GUI, but I want to do it programatically. I have tried the two following approaches.

Method #1:

raster = QgsProject.instance().mapLayersByName('itp')[0]
fcn = QgsColorRampShader()
fcn.setColorRampType(QgsColorRampShader.Discrete)
lst = [ QgsColorRampShader.ColorRampItem(22.6, QColor(0,255,0)),
QgsColorRampShader.ColorRampItem(23.2, QColor(255,255,0)),
QgsColorRampShader.ColorRampItem(100.0, QColor(255,0,0)) ]
fcn.setColorRampItemList(lst)
shader = QgsRasterShader()
shader.setRasterShaderFunction(fcn)
renderer = QgsSingleBandPseudoColorRenderer(raster.dataProvider(), 1, shader)
raster.setRenderer(renderer)
raster.triggerRepaint()

This method works, but the labels for the groups are empty (I do not know how to fix that), and I have to manually specify the expected ranges on the ColorRampShader, which is unconvenient. So I tried the following:

Method #2:

raster = QgsProject.instance().mapLayersByName('itp')[0]
colRamp = QgsGradientColorRamp (QColor(0,255,0), QColor(255,0,0), False)
renderer = QgsSingleBandPseudoColorRenderer(raster.dataProvider(), 1)
renderer.createShader(colRamp, QgsColorRampShader.Discrete, QgsColorRampShader.EqualInterval, 3)
raster.setRenderer(renderer)
raster.triggerRepaint()

The problem is this second method does not work at all...

...until I open the symbology properties and click accept in GUI. Calling triggerRepaint() does not help.

So, my question is:
Why Method #2 does not work? I could not find any funcion to force refresh the calculation inside the docs, and this method does not serve my purpose if I have to manually open every output file.

And, related to this: Is there any way to specify the number of decimal/significant digits of the group labels?


texts: Calculate the sine of a raster in the power of y / Python

texts:
I have this code and my aim is to calculate the sine of my raster in the power of 0.8.
import os
os.chdir('D:/NOA/Soil_Erosion/test_Project/Workspace/Input_Data_LS_Factor')
import rasterio
import math

data = rasterio.open('Slope_degrees_clipped.tif')
band = data.read(1) # array of float32 with size (3297,2537)
w = band.shape[0]
print(w)

h = band.shape[1]
print(h)

dtypes =data.dtypes[0]
band_calc = math.sin(band)**0.8 # the formula I would like to calculate

However, the following error pops up: only size-1 arrays can be converted to Python scalars.
May you know how I should fix this?
P.S1. I tried to vectorize it (np.vectorize()) but it does not work as it needs a real number.
P.S2. When I use the np.ndarray.flatten(band) the same error occurs.


texts: editor tracking updated in attribute table, but not in sql server database

texts:
When I enable editor tracking and make modification on an object, I can find the edit user and edit time, but I cannot find the data updated in the sql server database. Is it because of the versioning? If yes when will the data updated in the sql server database? If no, what is wrong with my sql server database?

The following pictures are attribute table in ArcMap and table value in sql server.

I am using ArcMap 10.2.2 and sql server 2008 r2.


texts:
ArcGIS Pro caches the contents of files even after the file has been replaced on disk.

How can this feature be disabled?